Содержание к диссертации
Введение
1. Численное моделирование пространственного магаитостатического поля. Программные комплексы КОМРОТ и KLONDIKE. 29
1.1 Введение 29
1.2 Дифференциальная постановка задачи магнитостатики с использованием векторного электрического потенциала 38
1.3 Задание векторного электрического потенциала Р и области его определения
1.4 Конечно-элементная аппроксимация 40
1.4.1 Метод Ритца 41
1.4.2. Метод Бубнова - Галеркина 44
1.5 Решение системы конечно-элементных уравнений 47
1.6 Численное моделирование распределения пондеромоторньгх сил 53
1.7 Программный комплекс КОМРОТ. Методические расчёты. 57
1.8 Интегральная постановка задачи. Программный комплекс KLONDIKE. Методические расчеты 65
1.9 Выводы 82
2. Численное моделирование пространственного квазистационарного электромагнитного поля в проводящих оболочках. Программный комплекс TYPHOON. 84
2.1 Введение 84
2.2 Постановка задачи расчёта квазистационарного электромагнитного поля с использованием векторного электрического потенциала
2.3 Принцип сведения многосвязного пространства к односвязному 97
2.4 Конечно-элементная аппроксимация 99
2.5 Вычисление интегралов по треугольным симплекс-элементам 105
2.6 Решение системы конечно-элементных уравнений 110
2.7 Программный комплекс TYPHOON. Методические расчёты. 117
2.8 Выводы 148
3. Математическая постановка задач синтеза магнитных систем 151
3.1 Введение 151
3.2 Основные понятия теории некорректных задач. Построение регуляризованного решения и выбор параметра регуляризации 161
3.3 Сведение задач синтеза магнитных систем к минимизации сглаживающего функционала
3.4 Численная реализация решения задач синтеза магнитных систем 166
3.5 Методические расчёты 169
3.5.1 Расчетная модель для анализа влияния ферромагнитных вставок на уровень гофрировки тороидального магнитного поля 171
3.5.2 Постановка задачи оптимизации влияния ферромагнитньк вставок
3.5.3 Численные результаты 175
3.6 Выводы 184
4. Результаты численного моделирования и формирования поля магнитных систем электрофизических устройств 186
4.1 Введение 186
4.2 Анализ поля ахроматичного поворотного магнита с азимутальной вариацией магнитного поля 197
4.3 Моделирование поля магнитов широкоапертурного спектрометра КОМБАС
4.4 Анализ и формирование поля септум-магнита 208
4.5 Формирование однородного поля спектрометрического магнита с постоянными магнитами 223
4.6 Синтез системы корректирующих катушек установки ИТЭР 222
4.6.1 Гармонический и статистический анализ ошибок поля электромагнитной системы установки ИТЭР
4.6.2 Синтез системы корректирующих катушек 230
4.7 Оптимизация магнитной системы экранирования инжекторов нейтральных пучков от поля рассеяния установки ИТЭР 246
4.7.1 Описание конструкции и критерии оптимизации системы магнитного экранирования основного инжектора нейтральных пучков 247
4.7.2 Метод и результаты численной оптимизации NB MFRS 249
4.7.3 Описание конструкции и критерии оптимизации системы магнитного экранирования диагностического инжектора нейтральных пучков 255
4.7.4 Результаты численной оптимизации DNB MFRS 256
4.8 Формирование изохронного магнитного поля в циклотроне DC 72. 273
4.9 Численное моделирование переходных электромагнитных процессов, вызванных срывами тока плазмы, в основных конструктивных элементах установки ИТЭР 292
4.9.1 Описание конструкции и расчётная модель вакуумной камеры 295
4.9.2 Описание конструкции и расчётная модель диверторной кассеты
4.9.3 Результаты расчёта электромагнитных нагрузок 298
5. Численная реконструкция магнитного поля в объёме по данным измерений компонент поля на замкнутой границе рассматриваемого объёма 307
5.1 Введение 307
5.2 Постановка задачи реконструкции магнитного поля 309
5.3 Алгоритмическая реализация. Комплекс программ MAGMAP 313
5.4 Результаты численных экспериментов 323
5.4.1 Анализ влияния ошибок на точность реконструкции магнитного поля
5.4.1.1 Восстановление поля по нормальной компоненте вектора индукции, измеренной на границе области (задача Неймана) 326
5.4.1.2 Восстановление поля по трем компонентам вектора индукции, измеренным на границе области (задача Дирихле) 332
5.4.2 Сопоставление результатов восстановления магнитного поля с данными магнитных измерений 334
5.5 Выводы 346
6. Определение положения и формы плазменного шнура в токамаках в режиме реального времени по данным внешних магнитных измерений 348
6.1 Введение 348
6.2 Постановка задачи реконструкции границы плазмы 352
6.3 Метод подвижных токовых нитей 356
6.4 Краткое описание электромагнитной системы установки GLOBUS -Ми системы магнитных измерений 364
6.5 Учет влияния вихревых токов 366
6.6 Результаты модельных расчётов. Численная реконструкция границы плазмы. 368
6.7 Выводы 393
Заключение 395
- Дифференциальная постановка задачи магнитостатики с использованием векторного электрического потенциала
- Интегральная постановка задачи. Программный комплекс KLONDIKE. Методические расчеты
- Конечно-элементная аппроксимация
- Основные понятия теории некорректных задач. Построение регуляризованного решения и выбор параметра регуляризации
Введение к работе
Современный уровень проектирования электрофизических установок и устройств предполагает широкое использование численного (вычислительного) эксперимента [1,2] с целью определения и оптимизации их основных параметров и характеристик. Во многих случаях в силу геометрической сложности установок и нелинейных свойств используемых материалов, численный эксперимент является практически единственно возможным источником получения необходимой информации. Он рассматривается как методологическая основа создания и изучения математических моделей исследуемого объекта с помощью вычислительных средств. Общие принципы его проведения наиболее полно сформулированы в работах [1,2].
Можно выделить следующие основные этапы вычислительного эксперимента: построение адекватной физическому явлению математической модели, выбор численных методов для её расчёта, разработка программного обеспечения, собственно расчёт и анализ результатов.
Характеризуя вычислительный эксперимент в целом, необходимо отметить его универсальность, что позволяет использовать эту технологию при исследовании различных объектов. Это обстоятельство характерно для математического моделирования, поскольку многие явления и процессы, протекающие в различных установках, имеют одни и те же математические модели. Вычислительный эксперимент значительно дешевле и доступнее по сравнению с натурным физическим экспериментом; его подготовка и проведение занимают меньше времени.
Ускорение темпов научно-технического прогресса, автоматизация научно-исследовательских и проектно-конструкторских работ с применением ЭВМ, а также связанное с ними сокращение сроков морального старения оборудования ставят задачу интенсификации процесса создания и освоения новой техники. При этом применение численных методов и соответствующего программного обеспечения даёт возможность использовать сложные нелинейные математические модели, охватывающие все существенные черты рассматриваемых физических процессов в широкой области изменения их параметров, проводить изучение устройств путём «проигрывания» их поведения в различных условиях, находить оптимальные параметры и режимы действующих или проектируемых установок.
В настоящее время проектирование и создание электрофизических установок в области ускорительной техники, физики высоких энергий и для исследований по проблеме управляемого термоядерного синтеза, которые относятся к сложным инженерно-техническим объектам, часто, с предельными параметрами, требует комплексных исследований, проведение которых невозможно без вычислительного эксперимента.
Одним из центральных вопросов численного моделирования физических процессов, определяющих критерии проектирования и условия эксплуатации электрофизической аппаратуры, является детальный расчёт пространственного распределения как стационарного, так и квазистационарного электромагнитного поля. Высокие стоимости крупных установок приводят к необходимости выполнения большого объёма расчётных работ, связанных с анализом возможных режимов эксплуатации и оптимизацией основных узлов установок. Расчётные работы продолжаются также на этапе экспериментальных исследований, обеспечивая интерпретацию результатов физического эксперимента и доводку устройств для получения заданных параметров и характеристик. Для устройств в области ускорительной техники, как правило, требуется проведение прецизионных расчётов.
Значительный вклад в создание базовых методов расчёта электромагнитных полей внесли российские учёные [9,10,12,13,15,24,26,30,36,38-40,49,52,131,305].
Создание в России на базе ведущих научных центров в области электрофизики исследовательских групп и коллективов, плодотворно взаимодействующих друг с другом, привело за последние десятилетия к весьма существенному прогрессу в численном моделировании физических процессов, протекающих в электрофизических установках. Во многих случаях в силу электромагнитного принципа, положенного в основу работы этих установок, результаты численного моделирования и анализа распределения электромагнитного поля являются определяющим фактором при их проектировании и создании.
В диссертационной работе решается комплекс задач, связанных с численным моделированием магнитного поля электрофизических устройств. Актуальность исследований определяется как предметной областью, так и необходимостью решения ряда методологических вопросов с целью обеспечения проведения численного эксперимента.
Диссертация является обобщением работ, выполненных в соответствии с научно-тематическими планами ФГУП «НИИЭФА им. Д.В. Ефремова».
Цель диссертационной работы состоит в комплексной разработке эффективных численных методов расчёта, формирования и реконструкции магнитного поля в электрофизических устройствах на основе единой методологической и алгоритмической базы. Решение этой задачи обеспечивает возможность проведения анализа, выбора параметров и оптимизации магнитных систем электрофизических устройств в процессе их разработки, проектирования и создания, что также является целью работы.
Первый этап вычислительного эксперимента связан с постановкой задачи и выбором математической модели. Естественно, что среди математически эквивалентных формулировок наибольшую практическую ценность представляют такие, которые приводят к алгоритмам, допускающим эффективную численную реализацию.
В настоящей работе для численного моделирования как стационарного, так и квазистационарного электромагнитного поля последовательно используется метод модифицированного скалярного потенциала [14,20,21,60-65]. Такая общность данного подхода, важная как с точки зрения методической, так и практической, позволяет его рассматривать в качестве единого подхода к численному моделированию электромагнитного поля в терминах скалярного магнитного потенциала. Это обстоятельство является принципиальным для разработки эффективных численных алгоритмов расчёта поля.
На основе данного метода, который относится к группе дифференциальных методов расчёта поля [10-21], сформулирована дифференциальная постановка задачи магнитостатики, ориентированная, в первую очередь, на детальное численное моделирование пространственных полей прецизионных магнитных систем. Практика численных расчётов [30,52,53] показывает необходимость наличия программного обеспечения, отвечающего различным подходам. В работе рассмотрена также постановка задачи магнитостатики на основе метода объёмных интегральных уравнений [118-120], который относится к группе интегральных методов [12,22-43].
Численное моделирование квазистационарного поля и вихревых токов в тонких проводящих оболочках с использованием метода модифицированного скалярного потенциала приводит к интегро-дифференциальному уравнению [126,127], в котором искомой является скалярная величина.
Математическая постановка задач синтеза магнитных систем, которые являются обратными по отношению к задачам анализа, базируется на использовании вариационного принципа регуляризации А.Н.Тихонова [213], естественного с точки зрения практических приложений. При этом задача синтеза магнитных систем может быть сведена к минимизации регуляризирующего функционала [250].
Численная реконструкция магнитного поля в объёме по данным измерений компонент поля на замкнутой границе этого объёма проводится на основе решения краевых задач [375-377] - внутренней задачи Дирихле и внутренней задачи Неймана.
Решение задачи магнитной диагностики, связанной с определением положения и формы плазменного шнура в токамаках в режиме реального времени по данным внешних магнитных измерений, строится на основе математической модели с использованием метода подвижных токовых нитей [429,430], обеспечивающего требуемую на практике точность определения границы плазмы и скорость вычислений. Данный подход приводит к интегральной постановке задачи для уравнения Грэда-Шафранова [321,322].
Основой численных методов решения многих классов уравнений является дискретизация задачи с последующим сведением, как правило, к системам алгебраических уравнений. Для дифференциальных и систем дифференциальных уравнений, обычно, дискретный аналог задачи строится на основе метода конечных элементов, который в настоящее время рассматривается как общий метод их численного решения [77].
В силу имеющих место с одной стороны ограничений на доступные вычислительные ресурсы (необходимую память и время решения задачи на
ЭВМ), а, с другой стороны, в силу значительного порядка (>106) решаемых систем, в общем случае, нелинейных алгебраических уравнений, являющихся дискретными аналогами уравнений поля, выбор метода их решения является важным моментом.
В результате проведённых исследований, как теоретического характера, касающихся вопроса сходимости двухслойных итерационных схем с несамосопряжённым оператором перехода [92], так и численных исследований [93,94,100] в качестве базового метода решения нелинейных систем алгебраических (конечно-элементных) уравнений выбран метод симметричной последовательной верхней релаксации (Symmetric Successive Overrelaxation Method - SSOR) [95] с процедурой полиномиального ускорения сходимости [86] на основе ВТ-процесса [96]. Такой итерационный алгоритм отличает достаточно высокая практическая эффективность, что, собственно, и предопределило успех в решении прикладных задач, требующих высокой степени их дискретизации.
Краевая задача Коши, к которой сводится задача расчёта вихревых токов в проводящих оболочках, применительно к расчёту электромагнитных процессов в токамаке-реакторе обладает характерными признаками жёсткости [172-179] и требует применения жёстко-устойчивых методов её решения. Эта задача решается на основе известной вычислительной процедуры С. Гира [176-178], осуществляющей автоматический контроль погрешности на шаге с изменением по этому критерию как степени метода, так и величины шага численного интегрирования, и на основе матричного разложения [168,172], позволяющего решить проблему ограничения на выбор величины шага численного интегрирования при его разгоне в жёстко-устойчивых методах.
Диссертация состоит из введения, шести глав и заключения.
В Главе 1 рассматриваются вопросы, связанные с разработкой алгоритмов численного моделирования пространственного магнитостатического поля электрофизических устройств.
В разделе 1.1 приведена общая постановка задачи магнитостатики, дан обзор современных методов расчёта магнитостатических полей, обращается внимание на необходимость наличия в практических расчётах программного обеспечения, отвечающего различным подходам. Обосновывается выбор для дифференциальных методов расчёта поля в качестве базового - метода модифицированного скалярного потенциала, допускающего эффективную численную реализацию.
В разделе 1.2 дана физическая интерпретация данного метода с введением понятия «магнитного листа» для реальных электрических контуров, образованных проводниками конечного сечения. В рамках этого подхода получен общий вид выражения для представления векторного электрического потенциала Р, который для однородного распределения плотности тока по сечению проводника определяется на классе кусочно-линейных функций.
В разделе 1.3 сформулированы общие принципы построения вектора Р и области его определения, подчинённые требованиям эффективной численной реализации. Предложенный способ задания Р и области его определения позволяет описывать на основе единой методологической и алгоритмической базы в качестве источников поля как катушки возбуждения различной геометрической формы, так и постоянные магниты.
Построению конечно-элементных аппроксимаций посвящен раздел 1.4. Здесь в рамках метода модифицированного скалярного потенциала рассмотрена эквивалентная вариационная формулировка задачи и отличающийся большей общностью проекционно-сеточный метод построения дискретного аналога на основе метода Бубнова-Галеркина.
В разделе 1.5 рассмотрен базовый итерационный алгоритм решения нелинейных систем алгебраических уравнений, полученных в результате дискретизации задачи, который строится на основе метода симметричной последовательной верхней релаксации с процедурой полиномиального ускорения сходимости с использованием ВТ-процесса. Описывается способ автоматизации выбора параметров итерационного алгоритма для управления нелинейным итерационным процессом на основе поведения векторной нормы поправки итерационных приближений.
В разделе 1.6 рассмотрен алгоритм численного моделирования распределений как поверхностных, так и объёмных пондеромоторных сил в нелинейных магнитных системах, согласованный с конечно-элементным представлением результатов расчёта поля.
Изложенный в разделах 1.1-1.6 подход составляет основу математического обеспечения комплекса программ КОМРОТ для расчёта пространственного магнитостатического поля электрофизических устройств, содержащих токовые, ферромагнитные как магнитомягкие, так и магнитотвёрдые элементы конструкций сложной геометрической формы с учётом нелинейных свойств используемых материалов. Краткое описание комплекса программ КОМРОТ приводится в разделе 1.7. Здесь также приведены результаты методических расчётов, имеющих самостоятельное практическое значение, связанных с определением электромагнитных параметров вентильного двигателя ВЭД-5, магнитная система которого является комбинированной (и охватывает собой общий случай), включая детальное численное моделирование магнитного поля одного из основных структурных элементов магнитной системы двигателя -«кейса», представляющего собой сборку постоянных магнитов в кожухе из медных и стальных пластин. Экспериментальные и расчётные данные максимального значения мгновенного электромагнитного момента двигателя совпадают в пределах погрешности 1% измерения этой величины. Интегральная погрешность совпадения экспериментальных и расчётных данных распределения поля вдоль поверхности «кейса» не превышает 10 %.
Интегральная постановка задачи магнитостатики относительно вектора намагничения на основе метода объёмных интегральных уравнений и её дискретизация с использованием метода коллокаций рассмотрена в разделе 1.8. Данный подход положен в основу комплекса программ KLONDIKE, краткое описание которого также приводится в этом разделе, включая методический расчёт, имеющий самостоятельное практическое значение и связанный с оптимизацией конструкции гексаполя на постоянных магнитах для ЭЦР-источника многозарядных ионов. Показана сходимость численного решения распределения поля гексаполя, имеющего секторную конструкцию, с ростом числа секторов к аналитическому решению для идеального гексаполя с непрерывным законом распределения намагниченности по кольцу.
В Главе 2 излагается численный алгоритм расчёта пространственного квазистационарного электромагнитного поля и вихревых токов в проводящих оболочках произвольной формы, расположенных в пространстве произвольным образом. Как показывает анализ, приближение тонких проводящих оболочек применимо к численному моделированию переходных электромагнитных процессов в установках типа токамак, ряд основных систем которых, такие как вакуумная камера, дивертор, тепловая защита, элементы силовых и опорных конструкций могут быть описаны с помощью такой модели, допускающей эффективную численную реализацию с использованием метода модифицированного скалярного потенциала.
В разделе 2.1 приведена общая постановка задачи расчёта квазистационарного электромагнитного поля в двухкомпонентной среде диэлектрик-проводник; рассмотрены условия практической применимости модели тонких проводящих оболочек.
Дифференциальная постановка задачи магнитостатики с использованием векторного электрического потенциала
Апгоритмическая реализация описания вектора Р и задание области его определения Q.P должны подчиняться ряду дополнительных требований, обеспечивающих эффективность проведения расчётов: - описание вектора Р должно быть удобным и охватывать все значимые с практической точки зрения случаи задания источников поля; - точность описания Р должна быть адекватна требуемой точности решения задачи; - важным является решение задачи минимизации объёмов QP, занятых источниками поля, с целью уменьшения трудоёмкости проведения расчётов. Заметим также, что в общем случае расчёта пространственного магнитного поля вектор Р всегда может быть построен так, что в выбранной системе координат имеет лишь две компоненты [9,69]. На базе изложенного выше подхода с использованием понятия «магнитных листов» [14] для реальных электрических контуров, образованных проводниками конечного сечения, был разработан единообразный способ описания (задания) вектора Р и области его определения [70,71]. Согласно этому способу, область определения Р представляется совокупностью многогранных элементов. В пределах каждого такого элемента вектор Р строится на основе линейной интерполяции, используя значения его компонент, определяемых в вершинах элемента согласно (1.2.8). В качестве примера, на рис. 1.3.1 - 1.3.2, наряду с областью QP, приведён ряд конфигураций обмоток (проводников), часто встречающихся в практических приложениях. Для кольцевых катушек, как частного случая плоских электрических контуров, указаны три способа задания области Qp, включающей в себя объём проводника: - область под катушкой (рис. 1.3.1 а); - область внутри катушки (рис. 1.3.1b); - область произвольной формы, опирающейся на поверхность проводника (рис. 1.3.1с). Первый способ требует дополнительно задания на плоскости антисимметрии z = 0 граничного условия в виде функциональной зависимости р = р(х, у). Второй способ нежелателен, если область QP пересекает область Qy. С другой стороны, ему отвечает минимальный объём области С1Р. Для конфигурации обмотки, соответствующей рис. 1.3.2а, область QP определяется направлением «проецирования» вектора Р вдоль указанных направлений П, и П2.
Наконец, в последнем варианте (рис. 1.3.2Ь) с учётом геометрических построений Р = Pxez в Dx и Р = -Р2ех в D2. В области перехода DX2 от Д к D2 вектор Р имеет обе компоненты Рх и Pz, которые могут быть определены в виде Рх = Р2 — и Ряд магнитных систем электрофизических устройств, в частности, применяемые в ускорительной технике магниты в качестве источников поля используют магнитотвердьіе материалы [72-74]. В первую очередь, это постоянные магниты, изготовленные из Sm-Co и Nd-Fe-B сплавов, кривая размагничивания которых В{Н) близка к прямой линии, что указывает на высокую степень постоянства вектора намагничения М [75]. Для описания магнитных свойств таких постоянных магнитов может быть принята модель «идеализированных ферромагнетиков» [3], согласно которой вектор намагничения М представляется в виде где вектор постоянного намагничения М0( ) является векторной функцией только пространственных координат, а величина магнитной восприимчивости к от напряжённости поля Н не зависит. В этом случае [76], используя тождественное выражению (1.1.3) представление где в безгистерезисных, магнитоизотропных средах и ju = ju0 (І + к) = ju0jur, //,. - относительная магнитная проницаемость, уравнение (1.2.1) преобразуется к виду Вектора Р и М0 могут быть объединены в единый вектор Д правой части уравнения (1.3.4). Такое объединение, впрочем, может быть выполнено неформальным путём, исходя из физической интерпретации вектора Д, указанной выше. Очевидно, что область определения Р при описании постоянных магнитов в качестве источника поля ограничена объёмом, занимаемом самими магнитами. В случае зависимости (1.3.3), имеющей общий характер, нелинейная связь между М и Н учитывается естественным образом, аналогичным учёту нелинейных свойств магнитомягких материалов (см. (1.2.1)). Таким образом, расчёт поля комбинированных магнитньгх систем, содержащих токовые, ферромагнитные магнитомягкие и магнитотвёрдые элементы конструкций, может быть выполнен на основе единой методологической и алгоритмической базы. Рассмотренный выше способ задания Р и области его определения был реализован в виде программных модулей, что позволило в значительной степени автоматизировать процедуру задания исходных данных, связанных с описанием источников поля. 1.4 Конечно-элементная аппроксимация В настоящее время метод конечных элементов является общим методом численного решения дифференциальных и систем дифференциальных уравнений [77,78]. Известно [77], что уравнения самого общего типа можно записать в слабой форме, допускающей обобщение от метода Ритца к методу Галеркина. В методе Ритца задача решения уравнения в операторном виде Lu = f, определённого на области Q, рассматривается как задача минимизации квадратичного функционала J (и) = (L и, и) - 2(/, и). Метод осреднения Галеркина является наиболее общим средством получения конечно-элементных моделей и представляет собой дискретизацию слабой формы {Lu- f,co)= \{Ьи - f)co dQ., редуцируя задачу к конечномерному случаю. С практической точки зрения важным является рациональный выбор весовых функций со((2), і = 1,2,... в соответствии с видом используемой конечно-элементной аппроксимации. Слабая форма дифференциального уравнения по сути является условием ортогональности невязки уравнения подпространству, порождаемому весовыми функциями соі, в соответствии с которыми происходит распределение невязки; идея поиска решения основывается на классической теореме о проекциях [78], фундаментальной для большинства методов осреднения.
Ввиду большей общности проекционно-сеточных методов [79] для численного решения прикладных задач, рассмотрение эквивалентной вариационной формулировки представляет, в основном, методологический интерес в рамках использования метода модифицированного скалярного потенциала Метод Ритца основан на том факте, что решение уравнения (1.2.1) совпадает с точкой минимума квадратичного функционала следующего вида где %\Н2)- величина, пропорциональная плотности энергии магнитного поля ju , для которой уравнение Эйлера для данной вариационной задачи имеет вид [80] которое, как можно убедиться, совпадает с (1.2.1). Здесь для определенности выбрана декартовая система координат х, у, z. На границе S = Si и S2 J S3 расчетной области V могут быть заданы граничные условия, следующие из условий симметрии антисимметрии магнитного поля или учета асимптотического поведения поля на больших расстояниях от локализованных внутри расчетной области источников поля [81,82]: где г{х, у, z) - радиус-вектор, проведенный из начала координат в точку наблюдения на границе S3 , q{x,y,z) - заданная функция пространственных координат. Эти краевые условия должны налагаться на класс функций, среди которых ищется минимум функционала (1.4.3). Вариационный подход позволяет исключить из специального рассмотрения некоторые граничные условия, если они являются естественными [77], т.е. автоматически удовлетворяются минимизирующей функцией. Предполагая зависимость величины магнитной проницаемости /л только от пространственньгх координат, что является оправданным при применении методов линеаризации для численного решения исходной краевой задачи, рассмотрим функционал вида При таком выборе функционала соответствующее ему уравнение Эйлера совпадает с (1.2.1), а естественное граничное условие имеет вид Условия (1.4.6) и (1.4.7) при совпадении единичных ортов ег =еп являются частными случаями (1.4.9), так что на минимизирующую функцию налагается только условие первого рода (1.4.5), которому легко удовлетворить на практике. Как отмечалось выше, можно построить вектор Р, локализованный вне области, занятой ферромагнетиком, так, что произведение juP обладает свойством а Рп = 0 на внешней границе расчетной области S3. Вьшолняя в (1.4.8) интегрирование по частям и переходя к поверхностным интегралам, можно получить Очевидным представляется соотношение.
Интегральная постановка задачи. Программный комплекс KLONDIKE. Методические расчеты
В процессе решения системы уравнений (1.8.6) значение jur корректируется каждый раз после получения приближения для М (метод линеаризации). Этот процесс продолжается до получения самосогласованного решения с заданной точностью. Для решения системы уравнений (1.8.6) используется метод простой итерации и ВТ-процесс ускорения сходимости итерационного процесса. Использование процедуры ускорения существенно сокращает время решения задачи. Данный подход, использующий метод объемных интегральных уравнений, положен в основу математического обеспечения комплекса вычислительных программ KLONDIKE [120, 121]. Основными особенностями этого комплекса являются: - использование точных аналитических решений для вычисления поля от однородно намагниченных элементов [119] и элементов с постоянным вектором плотности тока [122], ограниченных плоскими гранями; - использование эффективных численных процедур для вычисления поля от дугообразных токовых элементов с постоянным сечением [122]; - возможность одновременного учета различных магнитных материалов, в том числе, нелинейных. Комплекс состоит из набора независимых программных модулей, связанных единой файловой структурой, основными из которых являются: - препроцессор, позволяющий осуществлять ввод исходных данных в виде информационных файлов и визуализацию вводимых данных; - генератор системы алгебраических уравнений; - процессор, обеспечивающий решение нелинейной системы алгебраических уравнений; - постпроцессор, осуществляющий обработку и анализ результатов численного моделирования - нахождение компонент М, Н и В, вычисление пондеромоторных сил, действующих на элементы конструкций, формирование таблиц выходных данных. В последнем случае постпроцессор может быть использован, в том числе, для обработки результатов численного моделирования, полученных на основе комплекса программ КОМРОТ. В частности, для анализа распределения поля рассеяния магнитных систем. Комплекс программ KLONDIKE является эффективным инструментом для анализа «безжелезных» магнитных систем, в которых отсутствуют конструктивные элементы из ферромагнитных магнитомягких материалов. Верификация комплекса KLONDIKE проводилась в процессе численных исследований магнитного поля гексаполя [123] ЭЦР - источника многозарядных ионов с частотой накачки электромагнитной волны 14 ГГц [124].
В результате этих исследований была выполнена оптимизация массогабаритных параметров гексаполя при сохранении требований к величине и распределению магнитного поля. Как известно [125], для устойчивого удержания плазмы в ЭЦР - источнике многозарядных ионов используется магнитное поле конфигурации «минимум В», возрастающее во все стороны от области, занятой плазмой. Такая конфигурация поля в ионном источнике является результатом суперпозиции полей, создаваемых коаксиальными катушками и мультипольной структурой. Система катушек представляет собой простую зеркальную ловушку с пробочным отношением 1.5-7-2 и обеспечивает возрастание поля в аксиальном направлении. Обычно, в качестве мультипольной системы используется гексаполь, изготавливаемый на постоянных магнитах типа Sm-Co или Nd-Fe-B. Гексаполь создает возрастающее магнитное поле в радиальном направлении. Гексаполь представляет собой набор магнитных колец секторной конструкции, вектор намагниченности которых М лежит в плоскости, перпендикулярной продольной оси. Непрерывный закон распределения намагниченности по кольцу в азимутальном направлении М = М0 {ёу cos 3 p-er sin 3 р), (1.8.6) где er, ву - орты полярной системы координат, соответствующий случаю идеального гексаполя, в реальной конструкции аппроксимируется ломаной линией в соответствии с числом секторов, набираемых из однородно-текстурованного материала. Используемые в настоящее время материалы позволяют при описании их свойств с весьма хорошим приближением использовать модель постоянной намагниченности М0 = const во всем объёме постоянных магнитов. Для случая идеального гексаполя при заданных внутреннем радиусе Rx и внешнем радиусе R2 зависимость магнитной индукции от координат в представлении В-Н + АтиМ0, где Н - напряженность магнитного поля, имеет вид [123]: Заметим, что выражения (1.8.7) приведены в абсолютной (гауссовой) системе единиц, принятой в [123]. Анализ выражений (1.8.7) позволяет сделать следующие выводы: - в апертуре гексаполя {г Rx) магнитная индукция В г2; і— тек - в случае R2 j3Ri,B окрестности точек r = R{, ср- —, к = 0,1,2,...,5 рабочая точка на кривой размагничивания В = В(Н) лежит в третьем квадранте, что для современных магнитных материалов может рассматриваться как проявление нежелательных эффектов размагничивания внутренних слоев гексаполя его наружными слоями; - без учета эффектов размагничивания максимально возможное значение поля в апертуре гексаполя Вп = Втяк =1.5В0, где В0=АлМ0 - остаточная индукция, достигается на его внутренней поверхности при условии R2 » R{. В частном случае при R2 = J3RX - ВП=В0. Анализ краевых полей гексаполя требует проведения трехмерных расчётов. Оптимизация конструкции гексаполя при заданных длине и радиусе Rx, определяемых длиной и внешним радиусом ионизационной камеры, сводилась к варьированию радиусов R2 и числа секторов магнитных колец.
В качестве критерия рассматривалась эффективность использования материала постоянных магнитов (сплав Nd-Fe-B, В0 = 1.05 Тл) при сохранении требований к величине и распределению магнитного поля. В таблице 1.8.1 для соотношения R2=4bRx приведены результаты численных исследований зависимости от числа секторов величины магнитного поля Вп гексаполя, приведенной к максимально возможной величине В0. Как видно из таблицы 1.8.1, с ростом числа секторов полученное численное решение стремится к аналитическому решению для идеального гексаполя с непрерывным законом распределения намагниченности по кольцу. Оптимальным по технологическим соображениям представляется число секторов N = 24. Геометрические размеры гексаполя, состоящего из четырёх магнитных колец, набранных каждое из 24 секторов, приведены на рис. 1,8.1, Внешний диаметр крайних колец уменьшен, исходя из требуемой величины полного магнитного поля 0.5Тл на ЭЦР - поверхности нагрева плазмы с учётом поля системы катушек. Распределения магнитного поля в поперечных сечениях и в продольном направлении гексаполя в зависимости от радиуса приведено на рис. 1.8.2 - 1.8.3. На основе метода модифицированного скалярного магнитного потенциала разработана методика и алгоритм численного моделирования пространственных нелинейных краевых задач магнитостатики. В рамках этого подхода, используя понятие «магнитных листов» для реальных электрических контуров, образованных проводниками конечного сечения, и физическую интерпретацию векторного электрического потенциала как намагниченности этих листов, предложен единый способ описания источников поля - катушек возбуждения и постоянных магнитов. Таким образом, расчёт поля комбинированных магнитных систем, содержащих токовые, ферромагнитные магнитомягкие и магнитотвёрдые элементы конструкций, может быть выполнен на основе единой методологической и алгоритмической базы. На основе метода симметричной последовательной верхней релаксации и ВТ-процесса ускорения сходимости итерационных методов разработан эффективный итерационный алгоритм, отличающийся простотой реализации в сочетании с минимальными требованиями к вычислительным ресурсам. В результате на практике реализована возможность численного решения задач магнитостатики в областях сложной геометрической формы с большим количеством разнородных, нелинейных материалов. Таким задачам присуща высокая степень дискретизации (детализации), приводящая к плохо обусловленным системам, в общем случае, нелинейных алгебраических уравнений большого порядка ( 106).
Конечно-элементная аппроксимация
Простейшими среди всех конечных элементов являются так называемые симплексные [78, 164], которые связаны с использованием линейных интерполяционных функций в элементе. В двумерном евклидовом пространстве симплекс-элемент - это треугольник с тремя узлами, по одному в каждой вершине. Эти элементы используются в дальнейшем для построения конечно-элементной модели при численном решении дифференциального уравнения в частных производных (2.2.7), определённого на поверхности S. Также предполагается совпадение функций x N с базисными в пространстве решений при построении интерполяционных функций для приближённого представления искомой величины (нормальной компоненты вектора Р). Интерполяционные функции, которые также называются функциями формы, соответствующие двумерному симплексу, представляет собой полные полиномы первой степени вида ах + а2х + а3у. Разобьём поверхность S на совокупность конечного числа треугольных элементов 5,(е) с мерой пересечения, равной нулю: Введём локальную прямоугольную ортогональную систему координат x,y,z на каждом треугольнике так, чтобы ось z совпадала по направлению с нормалью п к плоскости треугольника, а оси х,у лежали в его плоскости. Тогда коэффициенты Ламе hx,hy,hz=l,n исходное интегро-дифференциальное уравнение (2.2.7) принимает вид: Рассмотрим первоначально конечно-элементный аналог дифференциального уравнения (2.4.1) на отдельном треугольнике, который затем с помощью достаточно простых отношений инцидентности может быть вложен в глобальную связанную конечно-элементную модель. Обозначим функции формы для каждого і -го узла треугольного элемента (е) через N . Заметим, что отношения инцидентности строятся на основе связи coi = и NJe . УМНОЖИМ уравнение (2.4.1) на JVJe) и проинтегрируем по поверхности элемента: Интегрируя по частям первый интеграл в (2.4.2), и применяя формулу Остроградского - Гаусса, получим: где lx - направляющий косинус нормали к контуру треугольника Z,(e) с направлением оси х. Преобразуя аналогичным образом второй интеграл в левой части уравнения (2.4.2), окончательно имеем: Примем, что функция Р на симплекс-элементе аппроксимируется с помощью тех же функций формы: где Рк -значение Р в к -ом узле треугольника; Nke) =—т -(а[е + b[e)x + c[e)y), при этом коэффициенты а[е\Ь1е\с[е определяются через значения глобальных координат вершин треугольника с учётом формул преобразования координат. Функции формы N[e) определены так, что они обладают следующим свойством: где 8кі (к,і = 1,2,3) - символ Кронекера, т.е. являются функциями с локальным носителем.
Для треугольного элемента наиболее распространенной является естественная система координат, определяемая тремя координатами Ц, L2 и Ьъ [78, 164]. Их значения дают относительные величины площадей треугольников, на которые разбит элемент L -координатами произвольной точки, принадлежащей конечному элементу. Преимуществом L -координат (координат площади [78]) является существование формул, упрощающих вычисление интегралов вдоль сторон Z/(e) элемента и по его площади (е) [165]: Выражение (2.4.5) позволяет осуществить связывание конечно-элементной модели. Первый интеграл в (2.4.4) преобразуется к виду: Здесь и далее индексам ink соответствует глобальная нумерация узлов; конечно-элементное уравнение строится по отношению к фиксированной глобальной узловой точке /. Применяя (2.4.7) к вычислению контурного интеграла в предыдущем выражении, получим: где 1$ представляет собой расстояние между двумя узлами і и к (длина соответствующей стороны треугольника). Второй интеграл в левой части (2.4.4) преобразуется к виду Возможно более точное вычисление этого интеграла с использованием, например, квадратурных формул Гаусса-Лежандра. Применяя (2.4.8) к вычислению поверхностного интеграла (2.4.12), получим: где 5ik - символ Кронекера. Второй интеграл в (2.4.11) с учётом соотношения (2.2.6), которое в принятой локальной системе координат x,y,z имеет вид J = V х Р =—ёх еу, преобразуется следующим образом: Выполняя дифференцирование подынтегрального выражения, представим двукратный интеграл в виде: окончательно имеем результат дискретизации уравнения (2.4.1) в форме: Подынтегральная функция в (2.4.15) сингулярна и применение квадратурных формул затруднительно. Для вычисления интегралов данного типа применялась методика, предложенная в работе [166], основанная на использовании понятия однородности функции.
В результате данные интегралы могут быть представлены в виде суммы интегралов меньшей кратности, которые в свою очередь, могут быть вычислены с использованием квадратурных формул для симплекс-элементов. В случае отсутствия вихревых и сторонних токов в начальный момент времени, J (г ,0) = 0, начальное состояние системы соответствует условию Начальное распределение Р(г,0) в случае резко выраженного скин-эффекта может быть найдено в результате решения системы алгебраических уравнений, исходя из начального условия где В[\0) - значения В \гк,0) в узловой точке к, а распределение поля внешних источников 5 0)(г,0) предполагается известным. 2.5 Вычисление интегралов по треугольным симплекс-элементам Возможны четыре случая взаимного расположения треугольников 5,(е) и (е} в пространстве при вычислении интегралов в (2.4.15): треугольники S и S полностью совпадают; треугольники 5(е) и (е} имеют одну и только одну общую сторону; треугольники S(e) и S(e имеют одну и только одну общую вершину; треугольники 5(е) и (е не имеют ни одной общей точки. В первых трёх случаях подынтегральная функция в (2.4.15) сингулярна. Рассмотрим вычисление интеграла в первом случае, когда треугольники 5,(е) и iS(e} совпадают. Выберем локальную декартову прямоугольную систему координат с началом в одной из вершин треугольника так, чтобы ось z совпадала по направлению с нормалью п к плоскости треугольника. Положение произвольной точки в пространстве относительно этой системы определяется координатами (x,y,z) = (x ,y ,zr), а ех=ех, ,еу=еу.
Основные понятия теории некорректных задач. Построение регуляризованного решения и выбор параметра регуляризации
Исторически развитие теории некорректно поставленных задач, возникающих в физике, технике и в других отраслях знаний, было связано с необходимостью создания теоретической основы для математической обработки и интерпретации результатов физического эксперимента [215]. Задачи синтеза по своей постановке существенно отличаются от задач интерпретации. Особенности математической постановки и выбор методов решения задач синтеза обусловлены последующим этапом проектирования и реализации систем и конструкций, т.е. связаны с необходимостью учёта условий физической и наилучшей конструктивной реализуемости результатов синтеза [261,262]. К первым относят условия, без выполнения которых невозможно реализовать решение задачи синтеза на практике. Таковым, например, для магнитных систем является требование гладкости распределения поля в области регулярности, для которой характеристики среды и характеристики поля конечны, непрерывны и имеют непрерывные производные. Ко второй группе относятся условия, учитывающие практическую реализуемость результатов синтеза в виде ограничений на сложность совокупности компонент системы. Условия физической и конструктивной реализуемости часто бывает трудно полностью разделить на практике. Оставаясь по своей исходной постановке некорректно поставленными обратными задачами математической физики, тем не менее, задачи синтеза требуют уточнения содержания и значимости отдельных составляющих понятия корректности. Во-первых, решение задачи синтеза следует понимать как возможность определения вектора совокупности компонент синтезируемой системы Z, для которого соответствующая ему характеристика u=Az (для магнитных систем, обычно, распределение поля) аппроксимирует заданную характеристику й с некоторой достаточной для практики точностью, т.е. M -M J. Следует иметь в виду, что в задачах синтеза жёсткие требования на точность аппроксимации заданной характеристики с учётом последующего этапа реализации результатов синтеза приводят к жёстким требованиям на технологические допуски изготовления и сборки элементов магнитной системы. Это вызывает как увеличение стоимости, так и ухудшение эксплуатационных свойств реализованных систем. Более того, часто в технических приложениях, требования на точность аппроксимации и степень выполнения условий конструктивной реализуемости приходится уточнять в самом процессе решения.
В задачах интерпретации при "точно измеренных" экспериментальных данных вопрос существования решения является априорно решённым положительно, так как измерения связаны с реально существующим физическим объектом. При математическом моделировании реального физического явления решение может не существовать, несмотря на априорное предположение о существовании точного решения, так как правая часть уравнения (3.1.1) получается на основании показаний физических приборов с погрешностью, определяемой принятой методикой измерений. Это одно из проявлений некорректности задач вычислительной диагностики. Во-вторых, вопрос единственности решения является принципиальным для задач вычислительной диагностики. Поэтому одной из главных математических проблем является выяснение условий, обеспечивающих единственность решения. В задачах синтеза возможны ситуации, когда целесообразно строить серии решений, а окончательный выбор решения задачи синтеза проводится с учётом условий наилучшей конструктивной реализуемости. Наконец, с точки зрения задач синтеза неустойчивость обратной задачи для уравнения (3.1.1) означает, что близкие характеристики системы могут быть реализованы на основе сильно различающихся по своим параметрам решений и, как следствие, по своим конструктивным особенностям, в том числе, не реализуемых или трудно реализуемых на практике. Таким образом, устойчивость решения задач синтеза также необходима и должна обеспечиваться за счёт включения в постановку задачи дополнительных условий, связанных с условиями наилучшей конструктивной реализуемости. Пусть вектор z - совокупность элементов синтезируемой магнитной системы, определяющая её заданные выходные параметры. В качестве варьируемых параметров (элементов) могут быть размеры, электрофизические свойства материалов, электрические параметры магнитных систем. Следует отметить, что в общей постановке задача проектирования оптимальных магнитных систем имеет комплексный характер и на практике типичной является ситуация, когда синтез магнитной системы сводится к оптимизации параметров выбранного типа конструкции. Обычно, вектор її - требуемое распределение магнитного поля в интересующей области, а А - оператор (матрица) прямой задачи. При заданном й вопрос о разрешимости задачи синтеза остаётся открытым, так как может не существовать элементов z є F, на которых Az = її, а точная нижняя грань функционала невязки Az - її \ может достигаться на элементах z, не реализуемых или трудно реализуемых на практике. Поэтому условие конструктивной реализуемости магнитных систем должно быть включено в математическую постановку задачи.
Во многих случаях такое условие может быть записано в виде квадратичного z\E d , где d - известное число. Тогда элемент z, на котором достигается точная нижняя грань функционала невязки, принадлежит подмножеству FjezF, Fj = {z : z 12Е d2} и удовлетворяет условию z=d2 [213]. Используя принцип взаимности для изопериметрических задач [263], эту задачу можно решать методом неопределённых множителей Лагранжа, т.е. искать элемент za, минимизирующий сглаживающий функционал АҐ [z, и$\ вида (3.2.3), а параметр а определять из условия гІ= 2. При этом стабилизирующий функционал \z\2E с практической точки зрения учитывает условия конструктивной реализуемости, а с математической - обеспечивает устойчивость решения задачи. Очевидно, что элемент za есть регуляризованное решение уравнения Az= її на множестве Fj. Три числовых параметра: 8 - величина погрешности; d — мера сложности элементов za и а-параметр регуляризации связаны двумя соотношениями: Задавая один из числовых параметров, два других однозначно определяются из (3.3.1). Таким образом, а может играть роль параметра, управляющего процессом синтеза: задавая для а некоторое множество значений и используя свойство строгой монотонности функций (f{a) и у/(а), определённых согласно (3.2.6), можно получить серию решений za, с различной степенью точности реализующих заданную функцию її. На основе такой серии решений оптимальное решение задачи синтеза может быть выбрано как наиболее приемлимое с точки зрения компромисса между требованиями достижения необходимой точности и выполнения условий конструктивной реализуемости. Следует заметить, что задача минимизации широкого класса функций большого числа переменных на множествах, определяемых ограничениями, может оказаться, как это часто бывает на практике, весьма сложной. Даже в наиболее простых задачах экстремум достигается в угловых точках границы множества условий, т. е. в точках, где нарушается дифференцируемость функции [258,264]. Кроме того, основные положения теории математического программирования и численные методы решения соответствующих экстремальных задач создавались в предположении, что задача корректна. С другой стороны, сглаживающий функционал М01 [z, и] можно ввести формально, не связывая его с вариационной задачей на условный экстремум, и строить регуляризирующий оператор путём решения задачи безусловной оптимизации. При этом а = а(8) согласно определению регуляризирующего оператора. Заметим также, что точная нижняя грань функционала невязки у2 = inf\Az-u\Ey являющаяся мерой несовместности уравнения Az = u, достигает своего значения на элементе z = Urn za, который по определению является квазирешением (см. 3.2.2). Задача синтеза может быть неразрешимой, если у превышает допустимый уровень 8тах, т. е. y Smax.
