На правах рукописи

БАТАРОНОВА Маргарита Игоревна

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ

МЕЗОСКОПИЧЕСКИХ СВЕРХПРОВОДЯЩИХ

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

КОНЕЧНО-ЭЛЕМЕНТНОГО АНАЛИЗА

Специальность 05.13.18 – Математическое моделирование,

численные методы и комплексы программ

АВТОРЕФЕРАТ

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

Воронеж – 2012

Работа выполнена в ФГБОУ ВПО «Воронежский государственный технический университет».

Научный руководитель Шунин Генадий Евгеньевич кандидат физико-математических наук, старший научный сотрудник, Воронежский государственный технический университет, профессор кафедры высшей математики и физико-математического моделирования

Официальные оппоненты: Ряжских Виктор Иванович доктор технических наук, профессор, Воронежский государственный университет инженерных технологий, зав. кафедрой высшей математики;

Пастернак Юрий Геннадьевич доктор технических наук, профессор, Воронежский государственный технический университет, профессор кафедры радиоэлектронных устройств и систем

Ведущая организация: ФГБОУ ВПО «Воронежский государственный университет»

Защита состоится 24 мая 2012 г. в 1100 в конференц-зале на заседании диссертационного совета Д 212.037.01 Воронежского государственного технического университета по адресу: 394026 Воронеж, Московский просп., 14.

С диссертацией можно ознакомиться в научной библиотеке Воронежского государственного технического университета.

Автореферат разослан 24 апреля 2012 г.

Ученый секретарь диссертационного совета Барабанов В.Ф.

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

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

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

МКЭ успешно применяется при решении задач в различных областях науки и техники. На его основе разработаны универсальные компьютерные системы инженерного анализа технических объектов, такие, как ANSYS, NISA, Cosmos/M, Maxwell, Elcut и др. Однако математические модели, используемые в этих системах, не включают уравнения электродинамики мезоскопических сверхпроводников, в силу чего их использование для моделирования МСЭМП невозможно. Требуется разработка дополнительных модулей к ним, либо применение специализированных систем конечно-элементного анализа. Для этих целей перспективно использовать наиболее мощную адаптируемую систему конечно-элементного мультифизического анализа COMSOL Multiphysics, а также разрабатываемый в Воронежском государственном техническом университете на кафедре высшей математики и физико-математического моделирования комплекс программ FEMPDESolver, предназначенный для численного решения скалярных дифференциальных уравнений в частных производных (ДУЧП) второго порядка методом конечных элементов c дополнительными условиями, учитывающими специфику электродинамики макроскопических сверхпроводников. С его помощью можно решать двух- и трехмерные задачи моделирования электрических и магнитных полей в многосвязных областях сложной формы в присутствии сверхпроводящих токонесущих элементов и рассчитывать электромеханические характеристики СЭМП. Для компьютерного моделирования МСЭМП требуется дальнейшее его развитие и расширение физико-математических моделей лежащих в его основе.

Данная диссертационная работа выполнена в рамках госбюджетных НИР:

Б5/07 «Моделирование топохимических и магнитомеханических процессов в многосвязных системах» (2007-2008 гг., № госрегистрации 01200707633), Б14/09 «Физико-математическое моделирование и исследование перспективных материалов, конструкций на основе титановых сплавов для авиационной и космической техники» (2009-2010 гг., № госрегистрации 01200952212), проводимых по заданию Федерального агентства по образованию в рамках тематического плана «Фундаментальные исследования», а также ГБ 2007.14 «Компьютерное моделирование криогенных магнитогравиинерциальных устройств», ГБ 2010.14 «Физико-математическое моделирование криогенных магнитогравиинерциальных устройств», и соответствует одному из основных научных направлений Воронежского государственного технического университета – «Наукоемкие технологии в машиностроении, авиастроении и ракетно-космической технике».

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

Для достижения этой цели были поставлены следующие задачи:

1. Рассмотреть основные типы СЭМП, методы их математического моделирования и возможности существующих универсальных систем мультифизического конечно-элементного анализа. Адаптировать комплекс программ FEMPDESolver для решения двухмерного уравнения Лондонов.

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

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

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

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

Тематика работы соответствует следующим пунктам паспорта специальности 05.13.18 «Математическое моделирование, численные методы и комплексы программ»: п. 2 «Развитие качественных и приближенных аналитических методов исследования математических моделей», п. 3 «Разработка, обоснование и тестирование эффективных вычислительных методов с применением современных компьютерных технологий» и п. 4 «Реализация эффективных численных методов и алгоритмов в виде комплексов проблемноориентированных программ для проведения вычислительного эксперимента».

Научная новизна работы. В диссертации получены следующие результаты, характеризующиеся научной новизной:

– дискретная математическая модель распределения плотности токов в системе токонесущих сверхпроводящих тел, построенная в рамках метода интегральных уравнений и учитывающая лондоновское проникновение магнитного поля в сверхпроводник ;

– метод регуляризации дискретной модели распределения плотности токов в системе сверхпроводящих тел, отличающийся сохранением хорошей обусловленности модели при увеличении числа тел;

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

– модифицированный комплекс программ FEMPDESolver 2.2, позволяющий численно решать ДУЧП эллиптического типа в многосвязных двухмерных и трёхмерных областях с заданными потоками градиентов потенциала и условием постоянства потенциала на внутренних границах, а также уравнение Лондонов в двухмерных областях;

– результаты вычислительных экспериментов, показавшие устойчивость подвеса мезоскопического сверхпроводящего кольца с постоянным магнитным потоком над (под) закреплённым мезоскопическим сверхпроводящим кольцом с постоянным магнитным потоком того же направления но другой величины;

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

Реализация и внедрение результатов работы. Процессор комплекса программ FEMPDESolver 2.2 зарегистрирован в ФГУП ВНТИЦ и внедрен в учебный процесс подготовки студентов специальностей «Техника и физика низких температур» и «Техническая физика» Воронежского государственного технического университета по дисциплине «Математические методы моделирования физических процессов».

Апробация работы. Результаты диссертационной работы докладывались и обсуждались: на II, IV-VIII Международных семинарах «Физикоматематическое моделирование систем» (Воронеж, 2005, 2007-2011); на V Международном семинаре «Компьютерное моделирование электромагнитных процессов в физических, химических и технических системах» (Воронеж, 2007); Воронежской зимней математической школе «Современные методы теории функций и смежные проблемы» (Воронеж, 2011); IV Международной научной конференции «Современные проблемы прикладной математики, теории управления и математического моделирования» (Воронеж, 2011); научнотехнических конференциях профессорско-преподавательского состава, аспирантов и студентов Воронежского государственного технического университета (2005–2011).

Публикации. Основные результаты диссертации опубликованы в 16 научных работах, в том числе 6 в изданиях, рекомендованных ВАК РФ. В работах, опубликованных в соавторстве и приведенных в автореферате, лично соискателю принадлежит: в [5,9,14,16] – модификация физико-математического, алгоритмического и программного обеспечения процессора комплекса программ FEMPDESolver; в [6] – аппроксимационные схемы и обобщение на трёхмерный случай; в [1-4,8,11-13] – проведение вычислительных экспериментов, участие в обсуждении результатов; в [15] – построение дискретных моделей и их численное исследование; [7,10,14]–анализ специализированных конечно-элементных программ.

Структура и объем работы. Диссертация состоит из введения, четырёх глав, заключения, списка литературы, включающего 185 наименований, и приложения. Основная часть работы изложена на 153 страницах и содержит рисунка и 25 таблиц.

ОСНОВНОЕ СОДЕРЖАНИЕ РАБОТЫ

Во введении обоснована актуальность темы, сформулированы цели и задачи исследования, определены научная новизна и практическая значимость результатов работы.

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

Конечно-элементный комплекс программ FEMPDESolver численно решает скалярные дифференциальные уравнения в частных производных эллиптического и параболического типов в двух- и трехмерных областях произвольной формы при заданных начальных данных, граничных условий 1, или 3 рода, скачка неизвестной функции и потока её градиента на разрезах, периодических условий. При разбиении области могут использоваться конечные элементы 1 и 2 порядка (треугольники, четырёхугольники, тетраэдры), в том числе изопараметрические, а также бесконечные элементы для решения задач открытого типа. Осуществляется визуализация результатов, получение интегральных и дифференциальных характеристик, графиков и таблиц. Развиваемая в диссертационной работе версия 2.2 этого комплекса отличается от предыдущих возможностью учитывать условия постоянства потенциала и заданный поток его градиента на замкнутых внутренних границах трёхмерных областей и решать уравнение Лондонов в двухмерных областях совместно с уравнением Пуассона. Для этого случая разработана специализированная пре/постпроцессорная оболочка со сторонним конечноэлементным разбивщиком Triangle. С помощью универсальной системы компьютерной математики Maple 14 разработан также программный модуль для решения системы интегральных уравнений Фредгольма 1-го и 2-го рода с логарифмической сингулярностью ядра.

Во второй главе в рамках метода интегральных уравнений на основе соотношения Лондонов для вектор-потенциала магнитного поля в сверхпроводнике формулируется математическая модель для распределения плотности сверхтока по объему лондоновских сверхпроводников и по поверхности мейснеровских сверхпроводников в виде системы интегральных уравнений Фредгольма 1-го и 2-го рода:

где k = 1..K, S0 – сечение тел, через которые протекает объемный сверхток с плотностью Jj, Сk – контур поперечного сечения k-тела, через который протекает поверхностный сверхток с плотностью Jk, K – число связных мейснеровских тел в системе, ядро K определяется формулой K = rr f ( m(r, z, r, z)), Фk – поток через контур k-тела, rk – радиальная цилиндрическая координата на поверхности этого тела, Aj – j-компонента вектор-потенциала внешнего магнитного поля. Ядро K в уравнении (1) имеет логарифмическую особенность при ( r, z ) = ( r, z), поэтому (1) относится к классу слабо сингулярных интегральных уравнений.

Регуляризация уравнений (1) физико-математической модели проводилась на основе метода выделения особенности и в результате была получена численная модель системы:

где D i, D j – шаги дискретизации по координатам, выбранным в сечениях тел, – матричные элементы ядра интегрального уравнения, а элементы K il, j получаются предельным переходом D i ® Dl, D j ® 0.

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

Исследование свойств модели для лондоновских сверхпроводников показало ее хорошую обусловленность при глубине l проникновения поля порядка размеров сечений тел, быстро ухудшающуюся с уменьшением l. Кроме того, ввиду сильной заполненности матрицы K решение модели даже при умеренном числе точек дискретизации требует значительных машинных ресурсов как по объему памяти, так и по времени счета. Практически применение модели ограничено двумя-тремя телами лондоновского типа.

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

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

По результатам моделирования вычислялись числа обусловленности nC, nE матрицы системы уравнений модели в чебышевской (С) и евклидовой (Е) нормах, сила Fш, действующая на пробный шар, а также силы, действующие на все тела системы. Кроме того, рассчитывались полные токи Ii на всех телах системы.

На основании полученных результатов по интегральным параметрам Fш, Ii методом Рунге-Ромберга определялась главная часть D ошибки модели и показатель р скорости сходимости.

Результаты исследования показывают линейную сходимость рассматриваемой модели по шагу дискретизации. При этом рассчитанные значения параметров в зависимости от числа шагов по экрану (шару) перестают меняться уже начиная с Nэ = 100, а параметры сходимости по Nк совпадают с полученными для всей системы. Это показывает, что для практических потребностей достаточно использования Nэ = 50100.

Зависимость числа обусловленности системы от числа точек дискретизации для рассмотренных систем приведена на рис. 1. Как видно, обусловленность модели резко ухудшается с ростом числа степеней свободы: nCµNo2.

Однако анализ результатов позволил обнаружить, что такой рост числа обусловленности связан с увеличением числа точек дискретизации по экрану (шару). Если выбрать данные при постоянном Nэ, то обусловленность системы перестает вообще зависеть от числа точек дискретизации (рис. 2).

Таким образом, значение Nэ = 50100 является не только достаточным для обеспечения хорошей точности модели, но и необходимым для обеспечения ее обусловленности. Для мезоскопической системы предлагаемый прием является фактически методом регуляризации модели.

Число обусловленности обнаруживает слабую корневую зависимость от числа тел в системе в регуляризованной модели (рис. 3). При этом относительная ошибка определения интегральных параметров оказывается ограниченной на уровне 0,1-0,3% (рис. 4).

Рис. 1. Число обусловленности Рис. 2. Число обусловленности модели от числа точек дискретиза- модели от числа точек дискретизации для системы различного числа ции по кольцам для системы различтел (числа у линий) ного числа тел (числа у линий) Рис. 3. Число обусловленности Рис. 4. Относительная ошибка регуляризованной модели от числа В результате регуляризованная модель допускает расчет с большим числом сверхпроводящих тел в системе.

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

В постановке ряда краевых задач, в частности, для уравнения Лапласа, встречается дополнительное условие постоянства значений функции вдоль некоторой границы Г1:

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

Данное условие допускает вариационную формулировку. Действительно, если рассматривать функционал в котором второе слагаемое равно произведению потенциала j и некоторой скалярной величины Q и отлично от нуля только на поверхности границы Г1, на множестве всех функций достаточной степени гладкости, заранее принимающих постоянное значение на Г1, то можно показать, что функция, доставляющая минимум этому функционалу, будет удовлетворять на W уравнению Лапласа 2u = 0, на границе, за исключением Г1 – естественному граничному условию u n = 0 (Неймана), и интегральному соотношению (Здесь нормаль n ' направлена от границы Г1 внутрь области W.) Из методов численного решения краевых задач, допускающих вариационную формулировку, наиболее эффективным является метод конечных элементов (МКЭ). Рассмотрим особенности применения МКЭ к решению краевой задачи для уравнения Лапласа с учетом дополнительных условий типа постоянства функции (2) и соотношения (3).

Условие u G = const приводит к тому, что все узловые переменные {ui}, соответствующие узлам на G1, заменяются одной ассоциированной переменной u*.

Новую систему алгебраических уравнений, учитывающую данные дополнительные условия (2) и (3), можно получить простой модификацией системы, где этого учета нет. Представим ключевые моменты этого алгоритма:

1) каждому отдельному интегральному условию вида (3) ставится в соответствие ассоциированная переменная u* – в реальности узловая переменная, связанная с одним из узлов на G1; в глобальной системе уравнений резервируется соответствующее ей уравнение, под исключенные степени свободы место не выделяется, 2) цикл по всем конечным элементам – формирование локальных систем происходит обычным образом – как в стандартном МКЭ, 3) если хотя бы один узел обрабатываемого конечного элемента принадлежит поверхности G1, данная локальная система подвергается перестройке по одной из представленных выше схем – в зависимости от того, сколько узлов этого элемента лежит на G1, 4) локальные системы добавляются в глобальную систему дискретных конечно-элементных уравнений по обычным правилам – в соответствии с локальной и глобальной нумерацией узлов, причем уравнения, которые соответствуют ассоциированной степени свободы u*, добавляются в соответствующее уравнение глобальной системы, 5) на завершающем этапе в правой части уравнения глобальной системы, которое соответствует u*, записывается соответствующее значение величины интеграла (3) со знаком минус (–Q).

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

В итоговой системе будет присутствовать уравнение, явно отвечающее за выполнение соотношения (3). Рассмотрев структуру вкладов в это уравнение, можно обнаружить помимо основных интегралов – по граням, целиком лежащим на G1, от u n, –дополнительный вклад в полный интеграл дают посторонние слагаемые – линейные интегралы по граням тетраэдров, имеющим только один или два узла на поверхности G1.

Если считать функцию u/n строго непрерывной в каждой точке сетки, то при суммировании вкладов соседних элементов, примыкающих к поверхности G1, посторонние интегралы будут взаимно уничтожаться. Поэтому фактически окончательно формировать уравнение глобальной системы, соответствующее ассоциированной степени свободы u*, будут поверхностные интегралы по граням, целиком принадлежащим поверхности G1.

Но непрерывность нормальной производной выполняется только для эрмитовых тетраэдров, которые имеют слишком большое число степеней свободы и поэтому редко используются на практике. Лагранжевы же тетраэдральные конечные элементы обеспечивают только непрерывность функции, но не частных производных. Поэтому градиент j является кусочно-непрерывной функцией, оставаясь непрерывной функцией внутри каждого элемента (постоянной функцией для тетраэдров 1-го порядка и линейной – для 2-го порядка).

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

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

В целом можно сделать вывод, что результаты аппроксимации интегрального условия (3) будут зависеть от качества конечно-элементной сетки, а также её регулярности. Отсутствие плохих (вытянутых) элементов и регулярная структура сетки обеспечивают плавное изменение u, а значит, более точное вычисление градиента. С другой стороны, с измельчением сетки в МКЭ точность вычисления решения задачи должна повышаться, а значит, и точность учета дополнительного условия. В этой связи заметим, что уменьшение размеров конечных элементов, а также повышения степени интерполяционного многочлена, действительно ведет к убыванию лишних слагаемых, вносящих ошибку.

В качестве тестовой задачи для уравнения Лондонов рассмотрен сверхпроводящий шар радиуса R ~ lL, помещенный во внешнее однородное магнитноеr поле напряженностью H0.Требуется найти распределение магнитного поля H внутри и в окрестности шара с учетом лондоновского проникновения. Выбрав цилиндрическую систему координат и поместив ее начало в центр шара, искомое распределение можно определить, решив уравнения где lL – так называемая лондоновская глубина проникновения,r y = 2prA – функция потока, A = A( r, z ) eq – векторный потенциал. Так как H = rot A, то компоненты напряженности магнитного поля выражаются через y так:

При решении этой задачи методом конечных элементов расчетная область не может уходить на бесконечность, поэтому она выбиралась в виде прямоугольника a z b, 0 r с достаточно больших размеров (рис. 5). Стороны этого прямоугольника были разнесены на 20R по переменной r и 10R – по z. Как показал вычислительный эксперимент, при таких размерах влиянием границ на поле шара можно пренебречь. Наиболее густая сетка использовалась внутри шара и в его непосредственной близости (примерно 35 элементов вдоль радиуса шара). Тип конечных элементов – треугольники Лагранжа 1-го порядка. Граничные условия следующие: на оси z функция y равна (из-за y = 2prA), на остальных границах ввиду удаления от шара поле должно быть приближенно равно внешнему, поэтому там ставятся условия Неймана:

Расчет произведен для случая lL = 0,3R. График изменения напряженности вдоль линии j = p 2 показан на рис. 6. Линии равных значений для функции y приведены на рис. 7, а. Для наглядности также на рис. 7, б изображен случай выталкивания поля из сверхпроводника при l L ® 0 (взято lL = 0,07R).

шара в однородном поле Рис. 7. Линии уровня поля y для случаев lL = 0,3R (а) и lL = 0,07R (б) Для рассматриваемой задачи известно точное решение, определенное для компонент поля Hr и Hj. Поэтому анализ погрешности конечноэлементного решения выполнен для производных функции y. Для лагранжевых элементов, которые дают непрерывное решение для исходной функции и разрывное для ее градиента, точность вычисления и сходимость производных величин обычно значительно ниже аналогичных показателей самой функции.

Тем не менее, непосредственное сравнение численного и аналитического решений для y r и y z показало достаточно хорошее соответствие. Во всех точках области относительная погрешность не превышает 1%, за исключением поверхности шара, где она равна 2%, а также в малой окрестности оси z, где она достигает 10%. Такая невысокая точность вблизи оси z объясняется наличием множителя 1/r в подынтегральных функциях элементов матрицы жесткости:

Тем самым обычные элементы 1-го порядка не в состоянии аппроксимировать столь быстрое изменение при r ® 0. Можно попытаться уменьшить ошибку вблизи z, если применить конечные элементы более высокого порядка либо разбиение повышенной плотности. Но особого смысла в этом нет, поскольку в остальной области обеспечивается удовлетворительная точность.

Более того, применение сглаживающих процедур, например, среднеквадратическое приближение исходной функции y полиномами 5-6-го порядка, позволяет на порядок уменьшить ошибку вычисления производных.

Отметим также, что достаточно малые значения параметра lL плохо влияют на обусловленность глобальной матрицы жесткости. Этому же способствует присутствие в подынтегральных выражениях множителя 1/r, что приводит к резкому увеличению (от нескольких до десятков раз) итераций метода сопряженных градиентов. Использование предобусловленности в виде {ri rj }, где ri – r-координата i-го узла конечно-элементной сетки, позволяет заметно улучшить сходимость данного метода.

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

Моделируемый объект представляет собой кольцо (радиус поперечного сечения a, осевой радиус b = 2,8a ), в котором течет сверхток I, в присутствии шара радиуса R = 1,5a, центр которого находится на оси симметрии кольца на расстоянии d от плоскости кольца. Геометрия области в цилиндрической системе координат показана на рис. 8, а. Очевидно, задача имеет осевую симметрию.

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

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

Расчет поля осуществлен с учетом лондоновского проникновения поля внутрь шара. В то же время предполагается, что из толщи кольца поле вытеснено полностью. Математически задача формулируется в виде уравнений Лондонов и Лапласа, задаваемых в соответствующих подобластях, для функции потока y = 2prA (A – единственная отличная от нуля z-координата векторного магнитного потенциала). Граничные условия на поверхности шара могут учитывать как фиксированный поток (условие 1-го рода y = F0), так и на на рис. 8, б в виде распределения линий уровня.

Рис. 8. Расчетная область и распределение функции потока системы кольцо + шар Проведен анализ подъемной силы, действующей на шар со стороны поля, в зависимости от его положения для трех значений лондоновской глубины проникновения lL = 0; 0,2a; 0,4a при постоянном магнитном потоке F0. Соответствующие кривые показаны на рис. 9 сплошными линиями. Максимум силы наблюдается при d » 0,8a. В режиме постоянного тока характер зависимости не меняется, только максимум смещается вправо и соответствует d » 1,0a. Одна из таких кривых, полученная при lL= 0,2a, отображена пунктирной линией.

Рис. 9. Зависимость силы от положения шара для разных значений lL.

F0 – максимальное значение силы, достигаемое при lL = Рассмотрим взаимодействие двух одинаковых сверхпроводящих колец с круговым сечением (радиусы колец 2,8 мкм, радиус сечения 1 мкм) в случае так называемого «замороженного магнитного потока», когда поток через отверстие каждого кольца сохраняет своё неизменное значение при взаимном перемещении колец. Тогда, как известно, в каждом кольце будет устанавливаться ток Ii, создающий такой собственный поток, чтобы в сумме с потоком от другого кольца дать фиксированную величину Fi. Распределение поле вблизи колец может быть найдено через скалярный магнитный потенциал u на основе вариационной формулировки. Требуемое распределение u, а также значения I1 и I2, доставляют минимум функционалу с учетом дополнительных условий где G1 и G2 – поверхности разреза, закрывающие отверстия колец, ( u+ - u- ) – скачок потенциала. Стационарное значение данного функционала соответствует уравнениям 2u = 0 (уравнение Лапласа в области W);

u n = 0 (граничные условия на границе G);

Для соосно расположенных колец, когда геометрия имеет осевую симметрию, задача носит двумерный характер. Были получены зависимости энергии и силы, действующей на одно из колец, в зависимости от расстояния между плоскостями колец для разных значений сохраняющихся потоков F1 и F2.

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

Математически рассматриваемая задача ставится во всем трехмерном пространстве, но для реальных численных расчетов область должна быть ограничена. Как правило, граница замыкания определяется эмпирически, т.е. таким образом, чтобы изменением энергии при увеличении размеров области можно было пренебречь. В данной серии расчетов внешняя граница области была выбрана в виде цилиндра, ось которого совпадает с осью колец. Радиус цилиндра – 17 мкм, его высота – 34 мкм. Конечно-элементные вычисления проводились с использованием лагранжевых тетраэдральных элементов 2-го порядка, имеющих 10 узлов. Полное число степеней свободы – около 640 тыс.

В случае, когда потоки F1 и F2 равны и одного знака (при этом F1 = F2 = = 20,710–15 Вб), силы взаимодействия между кольцами носят характер притяжения. Расхождение между результатами вычислений энергии и силы в трехмерном и осесимметричном случаях лежит в пределах 1.5–2%. Это объясняется тем, что размер расчетной области в трехмерном случае был выбран меньше, чем в двумерном (из-за непомерного роста количества элементов/узлов). Как показал вычислительный эксперимент, с уменьшением отличия соответственных размеров трехмерной и двумерной расчетных областей указанные кривые сближаются.

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

Наиболее интересным представляется случай, когда потоки в кольцах разные по абсолютному значению (но с одинаковым знаком). Пусть, в частности, F1 = 20,710–15 Вб, F2 = 2F1. Вид взаимодействия в данном случае меняется с изменением расстояния. На малых расстояниях друг от друга (примерно до радиуса колец) между кольцами действует сила отталкивания. На малых расстояниях друг от друга (примерно до 3/4 радиуса колец) между кольцами действует сила отталкивания, быстро нарастающая при их сближении. При удалении колец она сменяется силой притяжения, достигающей максимума (при d = 1,7 мкм), а затем медленно убывающей до нуля. Зависимость энергии от расстояния характеризуется наличием потенциальной ямы (рис. 10).

Чтобы проверить наличие устойчивого положения и по другим степеням свободы, проведены расчеты при боковых и угловых смещениях одного кольца относительно другого. Они выполнены для расстояния между плоскостями колец 0,75 мкм, что соответствует минимуму энергии при осевых смещениях. Соответствующие графики изменения энергии в зависимости от малых боковых смещений и поворотов показаны на рис. 11 (а, б).

2. 2. 2. 2. Рис. 10. Зависимости энергии (a) и силы (б) от расстояния между кольцами d (F2 = 2F1) 2. 2. 2. Рис. 11. Изменение энергии в зависимости от бокового смещения (а) и поворота (б) одного из колец (случай F2 = 2F1) В обоих случаях для энергии имеем параболическую кривую с минимумом в нулевой точке (несмещенное положение). Сила (момент сил) при таких малых возмущениях носит возвращающий характер и нарастает линейно. Таким образом, для системы двух колец с одинаково направленными и не равными по модулю сохраняющимися связанными магнитными потоками существует трехмерная потенциальная яма, в которой система обнаруживает устойчивость относительно возмущений по всем степеням свободы.

В целях обнаружения возможного устойчивого левитационного состояния было также проведено трехмерное компьютерное моделирование взаимодействия тонкого и толстого кольца. Радиус поперечного сечения первого (толстого) кольца r1 = 1 мкм, второго (тонкого) кольца r2 = 0,2 мкм, Радиусы колец одинаковы и, как и ранее, равны 2,8 мкм. Результаты вычислений показали, что данная система является неустойчивой по угловым возмущениям.

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

ОСНОВНЫЕ РЕЗУЛЬТАТЫ РАБОТЫ

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

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

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

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

Основные результаты диссертации опубликованы в следующих работах:

Публикации в изданиях, рекомендованных ВАК РФ 1. Компьютерное моделирование сверхпроводящих систем методом конечных элементов / М.И. Батаронова, C.А. Кострюков, В.В. Пешков, Г.Е.

Шунин // Вестник Воронежского государственного технического университета. – 2007. – Т. 3. – № 8. – С. 25-27.

2. Компьютерное моделирование двух сверхпроводящих колец с постоянными магнитными потоками / М.И. Батаронова, C.А. Кострюков, В.В.

Пешков, Г.Е. Шунин // Известия АН. Серия Физическая. – 2008. – Т. 72. – № 9. – С. 1271-1274.

3. Трехмерный анализ взаимодействия двух токонесущих сверхпроводящих колец / М.И. Батаронова, С.А. Кострюков, В.В. Пешков, Г.Е. Шунин // Системы управления и информационные технологии. – 2009. – № 3.2. – С.

212-214.

4. Батаронова М.И. Компьютерное моделирование левитации мезоскопической сверхпроводящей сферы в неоднородном магнитном поле / М.И.

Батаронова, С.А. Кострюков, Г.Е. Шунин // Системы управления и информационные технологии. – 2009. – № 4. – С. 52-53.

5. Конечно-элементный комплекс программ FEMPDESolver / C.А. Кострюков, В.В. Пешков, Г.Е. Шунин, М.И. Батаронова и др. // Системы управления и информационные технологии. – 2010. – № 4(42). – С. 52-57.

6. Учет условия постоянства неизвестной функции в конечноэлементном комплексе программ FEMPDESolver / М.И. Батаронова, С.А. Кострюков, В.В. Пешков, Г.Е. Шунин // Вестник Воронежского государственного технического университета. – 2010. – Т. 6. – № 11. – С. 227-230.

7. Шунин Г.Е. Свободно распространяемые программные средства для конечно-элементного анализа в сети Интернет / Г.Е. Шунин, Я.В. Сбитнев, М.И. Батаронова // Физико-математическое моделирование систем: материалы II междунар. семинара. – Воронеж: ВГТУ, 2005. – Ч. 3. – С. 18-25.

8. Батаронова М.И. Решение уравнения Лондонов с помощью пакета FEMPDESolver / М.И. Батаронова, С.А. Кострюков, Г.Е. Шунин // Компьютерное моделирование электромагнитных процессов в физических, химических и технических системах: материалы V Междунар. семинара. – Воронеж:

ГОУВПО «Воронежский государственный технический университет», 2007.

– С. 165-168.

9. Физико-математические модели комплекса программ FEMPDESolver / М.И. Батаронова, С.А. Кострюков, В.В. Пешков, Г.Е. Шунин // Физикоматематическое моделирование систем: материалы IV междунар. семинара. – Воронеж: ГОУВПО «Воронежский государственный технический университет», 2007. – Ч. 2. – С. 110-118.

10. Шунин Г.Е. Физико-математическое моделирование сверхпроводящих электромагнитных подвесов / Г.Е. Шунин, М.И. Батаронова // Физикоматематическое моделирование систем: материалы V междунар. семинара. – Воронеж: ВГТУ, 2008. – Ч. 3. – С. 3-13.

11. Батаронова М.И. Конечно-элементный анализ сверхпроводящих электромагнитных подвесов // М.И. Батаронова, Г.Е. Шунин // Современные методы теории функций и смежные проблемы: материалы Воронежской зимней математической школы. – Воронеж: Издательско-полиграфический центр Воронежского государственного университета, 2011. – С. 41-42.

12. Батаронова М.И. Структура системы компьютерного моделирования сверхпроводящих электромагнитных подвесов / М.И. Батаронова, Г.Е. Шунин // Современные проблемы прикладной математики, теории управления и математического моделирования: материалы IV Междунар. науч. конф. – Воронеж: Издательско-полиграфический центр Воронежского государственного университета, 2011. – С. 6-8.

13. Препроцессор комплекса программ FEMPDESolver / М.И. Батаронова, С.А. Кострюков, М.В. Мощёнский и др. // Физико-математическое моделирование систем: материалы VII междунар. семинара. – Воронеж: ФГБОУ ВПО «Воронежский государственный технический университет», 2011. – Ч.

3. – С. 104-108.

14. Шунин Г.Е. Разработка и моделирование сверхпроводящих электромагнитных подвесов / Г.Е. Шунин, М.И. Батаронова // Физикоматематическое моделирование систем: материалы VII междунар. семинара. – Воронеж: ФГБОУ ВПО «Воронежский государственный технический университет», 2011. – Ч. 3. – С. 109-128.

15. Батаронова М.И. Численный анализ интегральных уравнений описывающих распределение плотности токов в мезоскопических сверхпроводящих подвесах / М.И. Батаронова Г.Е. Шунин // Физико-математическое моделирование систем: материалы VIII междунар. семинара. – Воронеж: ФГБОУ ВПО «Воронежский государственный технический университет», 2012. – Ч.

3. – С. 24-50.

16. Программа. Процессор пакета программ FEMPDESolver 2.2 / М.И.

Батаронова, С. А. Кострюков, В.В. Пешков, Г. Е. Шунин. М.: ФГУП ВНТИЦ, 2010. № 50201001684.