автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.16, диссертация на тему:Математическое моделирование пространственных магнитостатических полей методом конечных элементов на векторной ЭВМ

Р Г Б ОД

'•; - П!'Т

ОБЪЕДИНЁННЫЙ ИНСТИТУТ ЯДЕРНЫХ ИССЛЕДОВАНИЙ

11-94-268 УДК 519.632.4+519.677

ЮЛДАШЕВА Марина Борисовна

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ПРОСТРАНСТВЕННЫХ МАГНИТОСТАТИЧЕСКИХ ПОЛЕЙ МЕТОДОМ КОНЕЧНЫХ ЭЛЕМЕНТОВ НА ВЕКТОРНОЙ ЭВМ

Специальность: 05.13.16 — применение вычислительной техники, математического моделирования и математических методов в научных исследованиях

Автореферат диссертации на соискание ученой степени кандидата физико-математических наук

Дубна 1994

Работа выполнена в Лаборатории сверхвысоких энергий Объединенного института ядерных исследовании

Научный руководитель-.

доктор физико-математических наук,

профессор

Е.Н.ЖИДКОВ

Официальные оппоненты:

доктор физико-математических наук,

профессор

11.НЛЗАИК1Ш

доктор технических наук

П.Л.ШЕЛАЕО

Ведущее научно-исследовательское учреждение: Российский Университет дружбы народе», г.Москва

Автореферат разослан

1994г.

Защита диссертации состоится 1994г. в {а

30

Лаборатории вычислительной техники и автоматизации Объединенного института ядерных исследований, г.Дубна.

С диссертацией можно ознакомиться в библиотеке ОШШ.

Ученый секретарь Совета / . 0

кандидат фиолко-математпчсскнх наук [УЛУи С З.М.ИВАНЧЕНКО

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

Актуальность

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

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

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

Работы, положенные в основу диссертации, выполнены в соответствии з проблемно-тематическим планом научно-исследовательских работ Объединенного института ядерных исследований.

Целью работы является создание компьютерной модели спектрометрического магнита СП-40, расчет пространственного магнитного поля СП-40 в направлении движения пучка частиц не только внутри магнита, но и в обла-:тн расположения регистрирующей аппаратуры и мишеней; разработка ал-'оритмов для численного моделирования магнитного поля во всем пространстве с использованием конечных и бесконечных элементов; разработка специальных алгоритмов решения нелинейных систем уравнений, возникающих зследствие дискретизации непрерывной задачи, для векторной ЭВМ; создание универсального комплекса программ для расчетов трехмерных магнитостати-1еских полей.

Научная новизна работы состоит в следующем.

1. Впервые осуществлено моделирование и анализ пространственного магнитного поля спектрометрического магнита СП-40 установки ЭКСЧАРМ во всем пространстве.

а) Проведено сравнение результатов расчетов с данными измерений внутри магнита.

б) Расчет магнитного поля в области, где расположена регистрирующая аппаратура и мишени, осуществлен по двум методикам; проведено сравнение этих методик.

2. Разработаны алгоритмы для численного моделирования магнитного поля во всем пространстве с использованием конечных и бесконечных элементов. Исследованы аппроксимирующие свойства бесконечных элементов.

3. Разработаны векторные алгоритмы решения системы нелинейных уравнений с разреженной матрицей, которая возникает как следствие конечно-элементной дискретизации непрерывной задачи.

4. В результате численных экспериментов установлено, что введение в итерационную схему решения системы нелинейных уравнений параметра гп€(0,1] позволяет существенно сократить число нелинейных итераций. Предложен алгоритм адаптивного выбора т„.

5. Создан универсальный комплекс программ, предназначенных! для численного моделирования пространственных магнитостатических полей, включающий разработанные алгоритмы, в том числе, для векторных ЭВМ.

Практическая ценность работы заключается в том, что результаты численного моделирования пространственного магнитного поля СП-40 нашли практическое применение при обработке результатов экспериментов на установке ЭКСЧАРМ. На основе предложенных в диссертации алгоритмов создан универсальный комплекс программ, который позволяет проводить расчеты магнитных систем с произвольной "разумной" геометрией во всем пространстве. Разработанные алгоритмы имеют самостоятельный интерес и могут быть использованы при численном решении задач экспериментальной физики и ускорительной техники, аналогичных рассматриваемым в диссертации.

Апробация работы

Основные результаты диссертационной работы докладывались на Международном совещании по программированию и математическим методам решения физических задач (Дубна, 1993 г.), на международной конференции по дифференциальным уравнениям EQUADIFF-8 (Братислава, 1993 г.), на Международном совещании International Workshop: Beam Dynamics and Optimization (1994 г., Санкт-Петербург), на XXIX и XXX конференциях Университета дружбы народов (1993 г.,1994 г.), на семинарах по вычислительной и прикладной математике ЛВТА ОИЯИ.

Публикации

По теме диссертации опубликовано 6 работ, в том числе в журнале "Математическое моделирование", в трудах совещаний и конференций, в сообщениях ОИЯИ.

Структура и объем работы

Диссертация состоит но введения, трех глав и заключения, содержит б таблиц ,2?рисунков, список литературы из 111 наименований и изложена на L13 страницах машинописного текста.

Во введении дано обоснование актуальности темы диссертации, форму-гируется задача магнитостатики во всем пространстве и кратко полагаются >сновные подходы к ее решению. Анализируются их достоинства и недо-:татки с точки зрения эффективности их использования при численном монтировании. В связи с этим отмечается, что необходимость нахождения рас-1ределения магнитного поля во всем пространстве обусловливает выбор дифференциального подхода для решения поставленной задачи: во-первых, он по-шоляет при незначительных затратах памяти ЭВМ выполнять расчеты поля : необходимой точностью внутри ферромагнетика; во-вторых, для более точ-гого учета условия на бесконечности используется техника бесконечных эле-яентов. В отличие от комбинированного подхода, использующего граничное гнтегральное уравнение, в этом случае сохраняется свойство симметрии матрицы, которая получается в результате дискретизации непрерывных ура-1 нений методом Галеркина.

Далее во введении приводятся цель, научная новизна, отмечается практи-геская ценность работы. Кратко рассмотрено содержание глав диссертации [ сформулированы основные положения, выносимые на защиту.

В первой главе сформулирована дифференциальная постановка задачи от-[осительно двух скалярных потенциалов, общего — ф и редуцированного -!>. Пусть О, ¡г — область, занимаемая ферромагнетиком. Выберем некоторую »бласть А, внутри которой находится Пр. Границу Г2 обозначим через Го, а 1а = П/Пр. Г — граница между и Пл- Тогда имеем следующую задачу:

СОДЕРЖАНИЕ ДИССЕРТАЦИИ

div(iJ,Vrp) = 0, x^ilp-,

(1)

div(V<j>) = 0, х&ПА;

условиями на границах:

¡¿(дф/дп) = дф/дп + n-Hs, г£Г;

ф = ф+фБ, хет;

(2)

<£«0, жег0,

Функция ¡1 = /x(|Vi/>|) известна, причем 1 < ¡1, < ¡1 < ¡л" где ц, и ц* — определенные константы. Вектор Hs вычисляется по закону Био - Савара:

(3)

где fls ~ область, занимаемая обмоткой, J - заданный вектор плотности тока, |х — т/1 - расстояние между точками х и у. Для нахождения потенциала ф3 решается уравнение Лапласа внутри области 0,р :

д^О) = о, х е nF, (4)

с краевыми условиями:

~ = -n-Hs, а:€Г, J n-Hsds = 0. (5)

г

Вектора Вх = цУф, Ях = Щ, В2 = Я2 = Уф + Н°, Д, Я,е(Х2(П))3, i = 1,2, удовлетворяют неравенствам [7]:

|АВ,|<2^|ДЯ,|, (С)

лД-дД>|дД|2,

где дД = Д!^ - Д<2), Affi = - г = 1,2,. Верхние индексы 1 и 2 означают, что используется два различных значения векторов.

Эти свойства обеспечивают на основе теории сильно-монотонных операторов существование и единственность обобщенного решения задачи (1)-(2) [5]. При использовании метода конечных элементов для дискретизации непрерывной задачи (1)-(2) свойства (6) обеспечивают также существование и сходимость конечно-элементных приближений к точному решению [8].

В §1.1 описывается метод конечных элементов для дискретизации непрерывной задачи (1)-(2),(4)-(5).

Будем предполагать, что область Q имеет непрерывные по Липшицу границы Г и Го. В качестве конечных элементов Qj выберем выпуклые шестигранники. Разобьем О. на конечные элементы так, чтобы Cl = Uf2j , = 0, когда i^j, предполагая также, что каждая грань элемента Пj является либо подмножеством Го, либо является гранью другого элемента, а граница Г составлена из граней выбранных элементов.

Введем в пространстве с прямоугольной декартовой системой координат относительно переменной £ = (£', £3) линейные функции формы N'm для куба [—1,1]3 в виде:

N^(0 = (1/8)(1 + <rrol^)(l + <W2)(1 + <W3), m = 1,.... 8,

коэффициенты <ттд выбираются равными 1 или -1 чтобы была равна в вершине с номером т и равна 0 в любой другой вершине. Учитывая, что уб [—1,1]3 переходит в элемент Qj при использовании преобразования:

** = Ь = 1,2,3,

т = 1

це хк]!п, к = 1,2,3 - координаты вершины элемента с номером т, можно олучить базисные функции на элементе - Nrn(^(x)),m = 1, ...8. Решение а дач (1)-(2), а также (4)-(5) будем искать в виде:

'/'(■г) = Е^ЫОД- хеПр;

1

ф(-г) = £ ФjNJ{.rl .тбйл; (7)

Ф5Н = хейг,

3

1С через гр},ф}, ф^ - обозначены нензвестные значения потенциалов в узлах ространственной сетки, полученной в результате дискретизации области одставляя выражения для г}> и ф через базисные функции (7) в обобщенную ■ормулировку задачи (1)-(2) и учитывая краевые условия, получим уравнение:

1£ ф, £ фJJv^,i■v^т](m=($)

О,, ^ег г,епА Пл

-jN.il ■ + £ ф? JvNi ■ УЛ^Ю. г а А

потенциал фБ находится из уравнения:

£ ф* I- У7У_, <1П = -1N,11- (9)

П г Г

Интегралы по областям П/.', , входящие в (8) - (9), вычисляются при омощи суммирования вкладов от отдельных элементов В пределах какого такого элемента /4 считается постоянной, ее аргумент вычисляется в ентре. Под центром элемента понимается точка, каждая координата которой вляется средним арифметическим координат вершин шестигранника. В §1.2 даются теоретические основы метода конечных элементов. Проводится проверка выполнения условий теоремы о существовании и цшственности обобщенного решения задачи (1)-(2).

Приводится теорема о сходимости конечно-элементных решений к непре-ывному обобщенному решению.

В §1.3 дается краткий обзор методов решения дискретпзованных уравне ний. Дискретизация непрерывной задачи (1)-(2) приводит к системе нелинсй ных уравнений с разреженной матрицей:

ВДн = /, м = р(у). (ю;

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

ДпУ"+1 ~ У" = ~(F,yn - /), п = 0: 1, ..., (и;

Тп

где y0,F0, т0 - заданы. В магнитостатике обычно используется два вида этогс процесса:

1. Вл = Fn, т„=1;

2. Вп = F'n, где Fn - матрица Якоби.

В атом параграфе приводится оценка скорости сходимости итерационногс процесса (11) при Вп = F„, вытекающая из общей теории [9]. Из этой оценкг, не следует, что значение тп = 1, которое обычно используется на практике является оптимальным. В диссертационной работе проведены численные эксперименты по влиянию тп на скорость сходимости итерационного процесса (11) [1]. Предложен алгоритм адаптивного выбора параметра т„, который меняется в зависимости от поведения нормы невязки. Это позволило в 1.5 - 2 раза сократить число итераций.

В этом же параграфе кратко описаны основные методы решения систем нелинейных уравнений.

Независимо от способа линеаризации исходной нелинейной системы на каждой итерации процесса (11) приходится решать систему линейных уравнений.

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

В §1.4 дана краткая информация о наиболее известных пакетах программ, предназначенных для моделирования пространственных магнитостатических полей и использующих метод конечных элементов. Отмечен опыт не очень удачной векторизации multigrid-алгоритма для пакета программ MAFIA [10].

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

В §2.1 описываются разработанные векторные алгоритмы для решения дискретизованной задачи [1],[2].

В 1989 году для наиболее известной программы TOSCA [11], предназначенной для расчетов трехмерных электромагнитных полей, были апробированы

векторные алгоритмы, которые использовались в методе неполного разложения Холесского с сопряженными градиентами. Авторамп работы [12] был предложен алгоритм умножения конечно-элементной матрицы на вектор, основанный на выделении в матрице "псевдодиагоналей". Таким образом, процедура умножения разреженной матрицы на вектор сводилась к операциям с векторами длиной < N, где N - размерность решаемой системы, однако при этом существенно использовались особенности параллельно-векторной ЭВМ CRAY Х-МР.

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

В результате предлагаемого многоцветного упорядочения матрица линеаризованной системы становится блочной, при этом каждый блок имеет только одну отличную от нуля диагональ. Естественная нумерация используется для стандартного метода неполного разложения Холесского, а многоцветное упорядочение - для алгоритмов умножения матрицы на вектор и решения системы с предобусловливателем. При таком комбинированном подходе удается, во-первых, сохранить свойства предобусловливателя, а, во-вторых, в значительной мере использовать возможности векторного компьютера, поскольку указанные алгоритмы при многоцветном упорядочении сводятся к операциям с векторами длиной < N/с, где N - размерность системы, а с - максимальное число отличных от нуля коэффициентов матрицы в одной строке. Процедуры умножения матрицы на вектор и решения треугольных систем оа счет использования векторов такой длины выполняются, например, на векторной ЭВМ CONVEX С120 в 2-3 раза быстрее, чем по стандартным программам, оттранслированным с векторной оптимизацией.

Далее рассматриваются особенности векторизации предобусловленного метода сопряженных градиентов. Приводятся три различных способа нумерации неизвестных:

- естественное упорядочение;

- упорядочение по обратному алгоритму Катхилла-Макки;

- многоцветное упорядочение.

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

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

Особенностью задачи, рассматриваемой в диссертации, является необходимость точного учета условия на бесконечности. В связи с этим, в §2.2 приводится техника использования бесконечных элементов [3]. Даются теоретические оценки аппроксимирующих свойств бесконечных элементов и погрешности дискретизации в этом случае.

В §2.3 приводится оптимизация алгоритмов вычисления потенциала ф5 на границе области, занимаемой ферромагнетиком [4]. На основе однозначного представления Нз = где Нз -поле от обмоток, вычисляемое по закону

Био-Савара (3), предлагается универсальный экономичный алгоритм вычисления потенциала ф5 на границе ферромагнетиков.

Третья глава посвящена расчету [5] пространственного магнитного поля спектрометрического магнита СП-40 установки ЭКСЧАРМ и описанию комплекса программ МЭРЕЗБ [б].

В §3.1 формулируется физическая постановка задачи.

Спектрометрический магнит СП-40 активно используется в экспериментальных установках для исследования распадных процессов заряженных частиц высоких энергий, в частности, в установке ЭКСЧАРМ, важнейшим элементом которой он является.

Магнит СП-40 представляет собой диполь Ш-образной формы с внешними размерами 450x323x305 см3, апертурой 100x48 см2, длиной ярма магнита 190 см п с рабочей величиной магнитного поля в центре - 0.75 Тл.

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

В §3.2 приводятся результаты численного моделирования. Покаоано, что в пределах длины магнита вычисленные значения поля и результаты измерений отличаются не более, чем на 2-10-3. Расчеты также показывают, что модули относительных компонент поля Вх/ В0 и В2/В0 на торцах магнита не превосходят 2,7 % и 4 % соответственно. В удаленной области расчеты сделаны по двум методикам, приводится сравнения этих методик.

В §3.3 дается краткое описание разработанного комплекса программ

МБГЕЗБ, предназначенного для моделирования пространственных магнито-статических полей. Приводятся особенности используемой дифференциальной постановки, основных алгоритмов, реализованных в комплексе, структура комплекса и его возможности. Приводится пример расчета тороидальной магнитной системы, для которой требовалось определить распределение магнитного поля внутрн ферромагнетика.

В заключении приводятся основные результаты диссертации:

1. При дискретизации непрерывной задачи магнитостатики методом конечных элементов (МКЭ) получается система нелинейных уравнений. Введенный в итерационный процесс решения этой системы параметр тп€(0,1] позволил существенно сократить число нелинейных итераций.

2. На основе однозначного представления — — , где Н$- поле от обмоток, вычисляемое по закону Био-Савара, предложен универсальный экономичный алгоритм определения потенциала на границе ферромагнетика.

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

4. Разработана новая схема представления конечно-элементной разреженной матрицы.

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

6. Выполнены расчеты трехмерного магнитостатического поля спектрометрического магнита СП-40 установки ЭКСЧАРМ методом конечных элементов и проведено сравнение с данными измерений. Показано, что в средней части магнита основная относительная компонента поля Ву/В0 по результатам расчетов и данным измерений отличается не более, чем на 210"3.

7. С помощью двух методик: МКЭ и МКЭ с бесконечными элементами-выполнены расчеты трехмерного магнитного поля СП-40 в области, где расположены регистрирующая аппаратура и мишени. Расчеты поля, выполненные по этим методикам, отличаются не более, чем на « 10~3.

8. Создан универсальный комплекс программ МБГЕЗБ, предназначенный для расчетов пространственных иагнитостатическпх полей, включающий разработанные векторные алгоритмы, благодаря чему он может особенно эффективно использоваться на векторных и параллельно-векторных компьютерах.

Работы, положенные в основу диссертации:

1. Е.П.Жидков, М.Б.Юлдашева, О.И.Юлдашев. Векторные алгоритмы для решения трехмерных нелинейных задач магнитостатики, ОИЯИ, Р11-94-161, Дубна, 1994. Математическое моделирование, т.6, N7, 1994.

2. Е.П.Жидков, М.Б.Юлдашева, О.И.Юлдашев. Решение трехмерных нелинейных уравнений магнитостатики на векторной ЭВМ. Тезисы докладов XXX конференции УДН. М.: РУДН, 1994.

3. E.P.Zhidkov, M.B.Yuldasheva, I.P.Yudin, O.I.Yuldashev. Computation of the spectrometer magnetic field in the detectors region. Preprint of the JINR, Ell-94-111, Dubna, 1994.

4. E.P.Zhidkov, O.I.Yuldashev, M.B.Yuldasheva. Optimization of some algorithms for computer simulation of the nonlinear 3D magnetostatic fields, Proc. Int. Conf., Dubna, 14-19 June, 1993, in: Programming and Mathematical Techniques in Physics, World Scientific, Singapore, 1994.

5. Е.П.Жидков, М.Б.Юлдашева, И.П.Юдин, О.И.Юлдашев. Математическое моделирование пространственного магнитного поля спектрометрического магнита СП-40, ОИЯИ, Р11-94-160, Дубна, 1994.

6. М.Б.Юлдашева, О.И.Юлдашев. Комплекс программ MSFE3D для расчетов пространственных магнитостатических полей. Версия 1.2, ОИЯИ, Р11-94-202, Дубна, 1994.

Цитируемая литература других авторов:

7. S.J.Polak, A.Wachters, A.de Beer. Proc. 1-st Compumag Conf., - Oxford, 1976, pp. 19-27.

8. Х.Гаевский, К.Грегер, К.Захариас. Нелинейные операторные уравнения и операторные дифференциальные уравнения.- М.: Мир, 1978.

9. А.А.Самарский, Е.С.Николаев. Методы решения сеточных уравнений. М., Наука, 1978.

10. F.Krawczyk, T.Weiland. Use of a multigrid solver in the MAFIA module S3 for electro-and magnetoststic problems. IEEE Trans. Mag, 1990, Vol.26, N2, p.747-750.

11. J.Simkin, C.W.Trowbridge. Three dimensional non-linear electromagnetic field computations using scalar potentials. Proc. IEE, 1980, vol. 127, Pt. B, N6, pp. 368-374.

12. N.J.Diserens, A.R.Mayhook. Experience in the use of vector processors for 3D static analysis. IEEE Trans. Mag., 1990, v. 26, N2, pp. 831-833.

Рукопись поступила в издательский отдел 19 июля 1994 года.