Проекционно-сеточные методы для решения нелинейных эллиптических задач с дифференциальными операторами векторного анализа тема автореферата и диссертации по математике, 01.01.07 ВАК РФ

004611739

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

11-2010-98 На правах рукописи

ЮЛДАШЕВ Олег Ирикевич

ПРОЕКЦИОННО-СЕТОЧНЫЕ МЕТОДЫ ДЛЯ РЕШЕНИЯ НЕЛИНЕЙНЫХ ЭЛЛИПТИЧЕСКИХ ЗАДАЧ С ДИФФЕРЕНЦИАЛЬНЫМИ ОПЕРАТОРАМИ ВЕКТОРНОГО АНАЛИЗА

Специальность: 01.01.07 — вычислительная математика

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

2 8 0КТ 2010

Дубна 2010

004611739

Работа выполнена в Лаборатории информационных технологий Объединённого института ядерных исследований

Научные консультанты:

заслуженный деятель науки РФ Жидков Евгений Петрович

доктор физико-математических наук Иванов Виктор Владимирович

Официальные оппоненты: доктор физико-математических наук,

профессор Вабищевич Пётр Николаевич

Институт математического моделирования РАН

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

профессор Гулин Алексей Владимирович

Московский государственный университет им. М.В.Ломоносова

доктор физико-математических наук Сычевский Сергей Евгеньевич НТЦ "Синтез", ФГУП "НИИЭФА им. Д.В.Ефремова"

Ведущая организация: Казанский (Приволжский) федеральный университет.

Защита диссертации состоится "29" октября 2010 г. в 14 час. на заседании диссертационного совета Д 720.001.04 при Объединённом институте ядерных исследований (Лаборатория информационных технологий) по адресу: 141980, г: Дубна, Московской области, ул. Жолио Кюри 6.

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

Автореферат разослан " ¿А." СЕНТЯБРЯ 2010 г.

Ученый секретарь ///лд - тг 1Г , „

- Иванченко Иосиф Моисеевич

диссертационного совета

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

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

Проекционно-сеточные методы получили широкое распространение в практических приложениях, благодаря их успешному применению в области теории упругости в начале прошлого века. В настоящее время этот класс методов является наиболее распространённым в инженерных расчётах и с успехом используется в самых различных областях теории и практики. В связи с этим следует упомянуть А.Н.Крылова, Л.В.Канторовича, В.И.Крылова, И.С.Березина, Н.П.Жидкова, Г.И.Марчука, A.A.Самарского, А.В.Гулина, Н.С.Бахвалова и других авторов основополагающих книг по численным методам.

Отметим наиболее общие достижения за последние два десятилетия в развитии проекционно-сеточных методов для решения эллиптических краевых задач: это различные алгоритмы решения дискретизованных нелинейных уравнений на последовательности сеток; р— и hp— версии метода Бубнова-Галёркина для непрерывных и разрывных базисных функций, а также для векторных рёберных и граневых (facet) базисных функций; кроме того, это смешанные методы, методы декомпозиции области, методы расчёта с использованием комбинированных систем (т.е. систем с дифференциальными и граничными интегральными уравнениями). Все эти подходы к численному решению задач в различных областях знаний интенсивно развиваются.

В то же время появился новый класс сложных задач, к которому относятся расчёты крупных электрофизических установок, например, для физических экспериментов, готовящихся на Большом Адронном Коллайдере (LHC) в Европейском Центре Ядерных Исследований (CERN, г. Женева, Швейцария). Одной из таких установок является детектор ALICE [1-4]. Одновременно в проектах российских физических экспериментов главным фактором становится минимизация затрат на создание установок, что стимулирует применение оригинальных конструктивных решений и, как следствие, интенсивное использование математического моделирования с развитием известных методов расчёта и созданием новых алгоритмов. При этом сложность моделирования электрофизических устройств часто связана с необходимостью учёта поведения векторных полей различной природы.

Можно выделить следующие принципиальные трудности, возникающие при решении подобных задач с помощью обычных подходов:

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

'А.Н.Тихонов, В.А.Арсенин. Методы решения некорректных задач. М., Наука, 1979, с.18

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

алгебраические системы со специфическими свойствами возникают при использовании некоторых алгоритмов для точного учёта убывании векторного поля на бесконечности и применении некоторых схем метода Галёркина с разрывными базисными функциями. Для сложных задач такого типа также нет достаточно эффективных итерационных методов;

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

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

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

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

Цель работы

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

- получение обобщённой формулировки, соответствующей задаче;

- исследование свойств оператора обобщённой формулировки в гильбертовых

пространствах;

- разработка алгоритма дискретизации обобщённого уравнения.

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

В первую очередь, это сохранение точности вычислений в следующих случаях:

- когда расчётная область включает несколько несвязных областей, в которых заданы

нелинейные уравнения;

- когда на границах областей, в которых поведение векторных полей задаётся разными

уравнениями, должны выполняться определённые условия;

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

порядка;

- при учёте условий убывания на бесконечности для задач во всём пространстве;

- при вычислении общей карты поля.

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

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

- "тёплых" и "сверхпроводящих" магнитов со сложной геометрией и сильноменяющимися

полями;

- магнитных систем с заданным поведением поля в определённой области (решение

обратных задач).

Научная новизна

Среди нелинейных задач эллиптического типа выделяется класс задач, заданных во всём пространстве. Трудность решения задач этого класса заключается в необходимости учёта условия убывания решения на бесконечности. Обычный подход к решению таких задач состоит в решении последовательности нелинейных задач для определения удалённой границы с асимптотическими условиями. В ранних работах автора для точного учёта поведения решения на бесконечности в качестве граничного условия было предложено использовать граничное интегральное уравнение. В результате получается комбинированная система, содержащая нелинейные дифференциальные уравнения в дивергентной или вихревой формах в ограниченной области и граничное условие в виде интегрального уравнения. Видимо, независимо от исследований автора через двенадцать лет в CERN была разработана программа, которая использовала аналогичный алгоритм включения интегрального уравнения в общую расчётную схему. Кроме того, эта программа применялась для расчёта главного дипольного магнита Большого Адронного Коллайдера. В диссертационной работе впервые на основе теории монотонных операторов исследованы комбинированные системы относительно специальных проекций векторов на пространства функций с интегрируемыми в квадрате градиентами и на пространства вектор-функций с интегрируемыми в квадрате роторами. Впервые получены комбинированные системы, использующие только один граничный оператор на произвольной границе вокруг расчётной области, включающей подобласти с нелинейными и линейными дифференциальными уравнениями. Решение задач с медленно убывающими на бесконечности искомыми функциями с помощью таких систем экономит общее время компьютерных вычислений. Одна из таких систем использовалась для моделирования дипольного магнита синхротрона СПИН ОИЯИ с помощью метода конечных элементов, и полученные результаты хорошо согласуются с расчётами по конечно-разностной программе.

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

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

Одним из свойств проекционно-сеточных схем, которые используются для решения сложных эллиптических задач, является локальная аппроксимация. Для её апостериорной оценки существует несколько различных методов, где допустимая величина локальной погрешности задаётся из опыта решения подобных задач. В диссертации для задач магнитостатики эту допустимую величину предлагается определять на основе данных интерполяции измеренного с высокой точностью магнитного поля. В качестве эталона используются измеренные магнитные поля для установок L3 (GERN) и ЭКСЧАРМ (ИФВЭ, г.Протвино). Предложенные в диссертации новые характеристики локальной аппроксимации позволяют не только получать расчётные магнитные поля высокого качества, но и оценивать качество полей, измеренных на определённой сетке.

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

поддержанного грантом РФФИ 97-02-16765. Во втором алгоритме строится компьютерная модель спектрометрического анализирующего магнита по заданной поворотной силе, размеру рабочей области и некоторым другим параметрам. Предложенные алгоритмы существенно сокращают общее время построения компьютерных моделей магнитов с заданными свойствами.

Практическая ценность

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

1) комбинированные формулировки с одним интегральным оператором;

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

3) алгоритм контроля точности на основе единых характеристик для вычисленных и измеренных магнитных полей;

4) итерационный процесс формирования однородного магнитного поля за счёт уточнения поверхности ферромагнетика;

5) алгоритм построения компьютерных моделей спектрометрических дипольных магнитов по заданной поворотной силе, размеру рабочей области и другим параметрам,

- апробировались при моделировании сложных магнитных систем и могут быть использованы для расчётов подобных систем.

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

Относительно приведённых в диссертации результатов расчётов конкретных магнитных систем можно отметить следующее:

1) компьютерные модели проектов магнитных систем для установки ALICE были разработаны на стадии формирования концепции дипольного магнита, проектировавшегося в Лаборатории высоких энергий ОИЯИ. Проведённые расчёты, в частности, использовались для обоснования окончательного варианта проекта сверхпроводящего дипольного магнита, который вошёл в Дополнение к техническому проекту эксперимента (CERN/LHCC 96-32), а также для обоснования первых проектов «тёплого» дипольного магнита, включая вариант с конусной конфигурацией магнитопровода, который является самым большим в мире. Для большинства расчётов проводилось сравнение с результатами, полученными по другим программам или по данным измерений магнитного поля. Во всех случаях было получено хорошее согласование.

2) задачи, решённые в рамках грантов РФФИ 97-02-16765 и 99-01-01103, позволили обосновать возможность создания спектрометрического магнита с поляризующими наконечниками для проекта эксперимента с поляризованной мишенью в ИТЭФ. Проект такой магнитной системы является уникальным, поскольку предполагает её работу как в режиме анализирующего дипольного магнита, так и в режиме поляризующего магнита.

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

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

Основные результаты диссертационной работы докладывались на Международных совещаниях по проблемам математического моделирования, программированию и математическим методам решения физических задач (Дубна, 1983, 1993), X, XI, XV Всесоюзных и Всероссийском совещании по ускорителям заряженных частиц (Дубна, 1986,1988, Протвино, 1996), II Республиканской конференции "Интегральные уравнения в прикладном моделировании" (Киев, 1986), Всесоюзных конференциях по вычислительной физике и математическому моделированию (Волгоград, 1988, 1989), рабочем совещании "Методы и программы расчёта магнитных полей" (Протвино, 1989), Международных конференциях по дифференциальным уравнениям EQUADIFF -7, 8 (Prague, 1989, Bratislava, 1993), Международной конференции по численному анализу ISNA'92 (Prague, 1992), Конференции по ускорителям заряженных частиц РАС-95 (Dallas, 1995), рабочих совещаниях коллабораций ALICE (CERN, 1996), EXCHARM (Дубна, 1995, 1997) и PANDA (FZ-Juelich, 2004), Конференции-семинаре "Математические модели сложных систем" (Тверь, ТвГУ, 1999), Всероссийских конференциях по проблемам физики, химии, математики, информатики и методики преподавания (Москва, РУДН, 1999, 2000, 2002, 2006), Международной конференции "Тихонов и современная математика" (Москва, МГУ, 2006), Международной конференции "Математическое моделирование и вычислительная физика" (Дубна, 2009), на семинаре кафедры вычислительных методов (ВМК, МГУ, 2009), на семинаре научно-исследовательского вычислительного отдела НТЦ «Синтез» (НИИЭФА, С.-Петербург, 2009), а также на Учёных советах ОИЯИ, семинарах по вычислительной и прикладной математике, вычислительной физике ЛИТ (JIBTA) ОИЯИ.

Публикации

Соискатель имеет 67 научных публикаций. Основные результаты, вошедшие в диссертацию, опубликованы в 28 печатных работах, в том числе в 9 публикациях в журналах, рекомендованных Высшей аттестационной комиссией ("Дифференциальные уравнения","Математическое моделирование", "Журнал вычислительной математики и математической физики", "Вестник Российского Университета Дружбы Народов", "Nuclear Instruments and Methods"), а также в 19 публикациях в журнале "Краткие сообщения ОИЯИ", трудах Всероссийских, Всесоюзных и Международных конференций, в сообщениях ОИЯИ и других научных изданиях.

Некоторые из работ автора, на основе которых написана диссертация, представлены в электронной библиотеке публикаций ОИЯИ http://wwwl.jinr.ru/ , в электронных библиотеках институтов КЕК (Япония) и CERN (г.Женева) 2.

Структура и объём диссертации

Диссертация состоит из введения, четырёх глав, заключения и трёх приложений. Она содержит 21 таблицу, 92 рисунка, список литературы из 214 наименований и изложена на 340 страницах, включая приложения. В каждом параграфе принята своя нумерация лемм и теорем.

2http://www.slac.stanford.edu/spires/hepnames/

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

Во введении приводится краткий обзор литературы по проекционно-сеточным методам решения нелинейных эллиптических задач, в том числе с помощью hp— конечно-элементных адаптивных алгоритмов. Как известно, разработка проекционно-сеточных методов восходит к работам W.Ritz(1908), И.Г.Бубнова(1913) и Б.Г.Галёркина(1915), причём метод Бубнова - Галёркина (взвешенных невязок) обладает очень большой общностью. Отметим также книги по проекционно-сеточным методам С.Г.Михлина, Х.Л.Смолицкого(1965), М.А.Красносельского, Г.М.Вайникко и др.(1969), P.Ciarlet (1978), Л.А.Оганесяна, Л.А.Руховца(1979), Г.И.Марчука, В.И.Агошкова(1981), В.В.Шайдурова(1989), Ch.Schwab(1998), O.C.Zienkiewicz, R.L.Taylor(2000), P.Monk(2003), P.Solin(2006). В настоящее время численные методы для решения подобных задач (в основном, линейных), в том числе с уравнениями Максвелла, интенсивно развиваются в связи с пионерскими работами R.Courant(1943), H.Whitney (1957), Г.И.Марчука(1961), Р.П.Федоренко(1961,1979), В.А.Кондратьева(1967), О. А. Ладыженской (1970), I.BabuSka(1973,1981), F.Brezzi(1974), В.И.Агошкова(1977), Е.А.Волкова(1979), J.C.Nedelec(1980), A.Bossavit(1987) и других. Среди современных исследований можно выделить как работы отдельных авторов, так и целых коллабораций, представляющих свои публикации не только в периодической печати, но и в WWW. Далее кратко раскрыто содержание глав диссертации, её аппробация и структура.

Первая глава посвящена проекционно-сеточным методам решения краевых задач для векторов с интегральными граничными условиями, точно учитывающими убывание решения на бесконечности. Как известно из работ Р.П.Федоренко (1961) и Н.С.Бахвалова (1966), особенностью численного решения линейных эллиптических задач в дифференциальном подходе является быстрая сходимость метода простой итерации для быстро осциллирующих компонент решения задачи. В случае же интегрального подхода, как показано, например, в работе Е.П.Жидкова, Б.Н.Хоромского, О.И.Юлдашева (1981), в методе простой итерации быстрее сходятся гладкие компоненты решения задачи. В связи с этим комбинированный подход, который объединяет дифференциальный и интегральный подходы и позволяет точно учитывать условие убывания решения на бесконечности, продолжает вызывать интерес. Среди первых публикаций, посвящённых разработке комбинированных методов, можно отметить работы O.C.Zienkiewicz, D.W.Kelly, P.Bettess (1977), C.Johnson, J.C.Nedelec (1980), Е.П.Жидкова с соавторами (1982,1992), S.J.Salon (1985), П.Г.Акишина (1991) и других. По-прежнему разрабатываются прикладные пакеты и программы на эту тему, что свидетельствует об актуальности излагаемых здесь результатов.

В §1.1 доказывается теорема 1.1.1 о компонентах непрерывно-дифференцируемого вектора. Как известно, каждый вектор можно разложить по направлению единичного вектора и ему перпендикулярному. В силу теоремы Гельмгольца единичный вектор может быть градиентом, ротором или их суммой. Поэтому при описании компонент вектора в новом разложении получаем различные системы уравнений для его определения.

Во-первых, из теоремы 1.1.1 вытекает абстрактная формулировка внутренней краевой задачи

V • (ai(x, й)й) = f(x,u), х ей]

V х (a2(i,i)i) = F(i,u), xei2; (1) a(x)n(aiu ■ n) + fi(x)n x (a2u x n) = gu x 6 ЗП, (2)

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

V- (Ai(x,u) х и) = f(x,u), х 6 П; V(A2(x,u)-u) + Vx(A3(x,u)xu) = F(x,u), х € П; (3)

a(x)n((Ä1 х и) ■ п) + (Ä2 ■ й)) + ß{x)n х (Ä3 х u) х n) = g2, x 6 dü. (4)

Также из теоремы 1.1.1 следует формулировка задачи

V • (ai(x,u)u) + V • (^(х,«) х й) = f(x,u), х е П; V х (а2(х,й)и) + V{Ä2(x,u) ■ й) + V х (Ä3{x,u) х и) = F(x,ü), iSfi; (5) а(х)й((а]И • п) + {Äi хй)-п) + (Л2 • u)) + ß(х)п х ((а2и + Л3 х й) х п) = д3, х 6 ÖS1 (6)

Все три абстрактные задачи (1),(2), (3),(4) и (5),(6) при определённых коэффициентах и правых частях имеют смысл, поскольку уравнения каждой из них описывают поведение компонент непрерывно-дифференцируемого вектора.

Сформулируем абстрактные задачи для описания стационарных векторных полей во всём пространстве. Если дифференциальные уравнения (1) дополнить условиями

limM = 0, ¿=1,2, (7)

то имеет смысл задача (1),(7). А при условии на бесконечности

Шп fc^l = 0, ¿ = 1,3, lim 1^=0, (8)

|xj—.оо |x| |x|-.oo |z|

имеет смысл задача (3),(8). Аналогично, имеет смысл задача во всем пространстве (5),(7),(8). Для краткости дальнейшего изложения рассматривается частный случай задачи (1),(7) в следующем виде:

V • (ai{x,u)u) = f{x), V х (a2(x,u)u) = F(x), ie!l;

V-la^x)«) = /i(x), Vx(o1(x)u) = F1(x), ieni = (R3\n); (9)

n • (aiu — äiu) = 0, n x (a2S — äi«) = 0, хбГн öfi;

lim äiff = 0. |i|->00

Здесь граница Г, разделяющая области П и Пь является липшицевой, / 6 Ь2{П), /i € L2(Ch), F е (Ь2(П))3, Fi е (Ь2{Hl))3, все коэффициенты уравнений предполагаются непрерывно-дифференцируемыми в U U fii.

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

Частным случаем задачи (9) являются задачи магнитостатики. Здесь следует отметить, что в автореферате кандидатской диссертации автора 3 для двумерных задач

3 О.И.Юлдашев. Граничные интегральные уравнения в задачах магнитостатики. Автореферат диссертации на соискание ученой степени кандидата физико-математических наук. ОИЯИ, 11-87-111, Дубна, 1987

магнитостатики было предложено рассматривать граничные интегральные уравнения как специальные граничные условия (лемма 3), и на основании лемм 4,5 были обоснованы два метода включения граничных интегральных уравнений для векторного потенциала в конечно-разностную сеточную схему решения дифференциальной нелинейной задачи магнитостатики относительно векторного потенциала. По всей видимости, независимо от этих исследований, один из упомянутых методов включения, только в конечно-элементную схему, был реализован в виде специальной программы для решения подобного класса задач i. Впоследствии она использовалась для расчётов магнитного поля главного дипольного магнита Большого Адронного Коллайдера 5.

Следует отметить, что излагаемый в этой главе подход (§ 1.4) является более общим, чем методы, изложенные в работах M.J.Friedman(1984) и J.E.Pasciak(1984), а также более экономичным по затратам памяти ЭВМ, чем подходы из 5'6.

Для исследования системы (9) удобно ввести обозначения

^. = |Q1U, i6il, f=iQ2"' хеП> U = V-Z 1 а\й, х 6 fij, 1 äiü, x 6 Пь

Для i6i! будем использовать следующие зависимости: Z = i/(x, V)V и V = к(х, Z)Z.

При / = 0 для решения системы (9) рассматривается интегральное уравнение относительно векторного потенциала А, который вводится по формуле V = V х А.

Векторный потенциал удовлетворяет уравнению

RÄ=Ä~~ [ö(Ä)x V-dy = i,(x), 2 еП. (10)

47Г J Г

a

Здесь U(A) = (1 — v(x,V x y!))V x А, а г = — y\ - расстояние между точками x и у, А/(х) - заданный вектор, Aj 6 (L2(Q))3. Предполагается, что 0 < v, < v < 1, V, = const.

В § 1.2 исследуется интегральное уравнение (10) для векторного потенциала в гильбертовом пространстве

H(rottQ) = {и-, и е (Ь2(П))3, V х и € (Ь2(П))3} со скалярным произведением

(и, v)a = q(u, v) + (V х й, V х v),

где а = const, а > 0, ü, v е H(rot,il), (•, •) - скалярное произведение в (L2(П))3.

Через || • || и || • ||а обозначим соответствующие нормы в пространствах (L2{tt))3 и H(rot, П).

Основным условием, которому должны удовлетворять коэффициенты системы (9), является следующее неравенство:

\Ui~ü2\<{i-v.)\vx-v2\, х ей. (11)

Свойства оператора R : H(rot, П) —» H(rot, Г2) устанавливаются в следующей теореме.

4 S.Kurz, J.Fetzer, W.M.Rucker. Coupled BEM-FEM Methods for 3D field calculations with iron saturation. Proceedings of the First International ROXIE users meeting and workshop, 1999. CERN, open-2000-151, 09/Aug/2000.

5 S.Kurz, S.Russenschuck, N.Siegel. Accurate calculation of fringe fields in the LHC main

dipole. CERN, LHC Project Report 357,1999.

Теорема 1.2.1

Если выполняется условие (11), то

1) R - сильно монотонный оператор, то есть для Ai, А2 € H(rot,£l)

(Е(Аг) - R(Ai),A! - А2)а > XiWA - A2HI где a,x 1 - положительные константы;

2) R удовлетворяет условию Липшица, то есть для А1гА2 £ H{rot,U)

ЩА,) - R(A2)\\a < h\\A, - А2\\а,

где а,1\ - положительные константы.

В теореме 1.2.2 устанавливаются свойства производной Фреше оператора R:

R'(A)A+ = А+- QU',

где U's (Е- Z'{V))V+, Z'(V) = {dZi/dVj}iJ=1,2,3, V+ = V x A+. Формулируются условия при которых решение уравнения (10) можно получить с помощью итерационного процесса Ньютона.

Для дискретизации интегрального уравнения используется метод Галёркина. Вводится пространство

Jl{Q) = {и:йе (Wj1 (П))3, V ■ и = 0, и • й|г = 0},

где и • п\Т - след на границе Г = ЗП, а Wj(fi) - пространство Соболева со скалярным произведением

(u,v)ji = (V х и, V х v)

и норма ||m||ji = ||V х гГ||, где u,ve J:(Q).

В леммах 1.2.1 и 1.2.2 доказываются вспомогательные результаты об эквивалентности норм и свойств оператора Vx для элементов пространства Далее рассматривается

уравнение

PAR(A) - Af) = 0,

и доказывается теорема 1.2.3. Теорема 1.2.3

Если выполняется неравенство (11), то

1) Pj\R\ Jl{П) —> J^fi) - сильно монотонный и липшиц-непрерывный оператор, то есть для А\, Ai € J1 (П) выполняется неравенство

(Pj'(R(Ai) - R(A2)),A, - A2)j, > и.ЦА, - А2\\*г,

а также

- Я(Л2))|Ь. < (■2-к)\\А1 -X2\\ji.

2) В предположениях теоремы 1.2.2 оператор PjiR'(A) : J1 (fi) —> J1 (О) непрерывно обратим и удовлетворяет условию Липшица по А. При этом для А,А+ е J!(f!) выполняются следующие неравенства:

||Рл#(Л)Л+||л > „.||А+у,

||Pj.(#(A,) - Н&тл < - Л2|Ьь

где 7 = const.

Согласно методу Галёркина, уравнению (10) соответствует система

P^(fi(Am) - Aj) = 0,

где J^ - m-мерное линейное подпространство в Jl{Q) с нормой, индуцированной из ./'(П), базисом щ,..., йт, и

т

Ат =

«=1

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

Из определения нормы J^fi) следует, что уравнение (10) эквивалентно удобной для численной реализации комбинированной m-мерной системе

J i/V х Ат ■ V х Uidx + J п ■ V х Ui J ■ n)dSydSx = J{n x щ) • V x A/dS,

n г г г

i = l,2,..., m;

здесь UJm = (1 — Pjmz/)V x Am, а проекция -PjmfV x Am находится из специального уравнения. Jm - конечномерное подпространство в пространстве J(fi) = [и : й £ (^(П))3, V • й = 0}.

В § 1.3 исследуется интегральное уравнение для скалярного потенциала. В случае, когда F = 0, в задаче (9) можно ввести скалярный потенциал </> по формуле Z = V(p. Он удовлетворяет уравнению

T{<p) = <p(x)+-^Ju(<p)-V^dy = <pf(x), хеп. (12)

п

Здесь Ufa) = (k(x,Vip) — 1)V((3, a ipj - заданная функция, <р/ € Ь2(П). Предполагается, что 1 < к < к*, к* = const.

Пусть Gi(fi) - гильбертово пространство: Gi(ii) = {й : и = е Wjt^)}. через

: Gi(fi) —► Gi(fi) обозначим интегральный оператор

(А0и)(х) = V f u-V~dy.

J 4ят а

Известно, что оператор А0 самосопряжён и положительно определён. Кроме того,

Ai < Moll < < 1, Хк = const, At > 0, к = 1,2. Учитывая эти свойства оператора А0, в Gi(fi) удобно ввести скалярное произведение

(и, гГ), = (A^u,v). Для соответствующей нормы имеем

iMi ^ и, < -^г1 н "II ■

и

В пространстве Wj (П) введём скалярное произведение

(р, 1/0,5 = /?(</>,'(iO + iv^vvO,,

где /3 = const, (3 > 0, <р,ф £ И^П), с соответствующей нормой || • = (-, Далее предположим, что выполняются неравенства

\\Ui - U2\\ < {к' - 1)\\ZI-22\\, (13)

(U1-U2,Z1-Z2)>0, (14)

где Ut, Zi е (12(П))3, г = 1, 2.

Справедлива следующая теорема. Теорема 1.3.1 Если выполняются неравенства (13) и (Ц), то

1) Т является сильно монотонным оператором, то есть для любых <pi,>p2 G Wj(f2)

(T(tpi) - T(tp2),tp! - <p2)0 > xt\\tfi - vail % где P, X2 - положительные константы;

2) T является липшиц-непрерывным оператором, то есть для любых <fi,tp2 6 Wj(n)

где 0,12 - положительные константы.

Для решения уравнения (12) методом Ньютона определяется производная Фреше

T(v)<p+ = у>+ + ¿ J(V'(Z) - 1 )2+-4±dy,

п

где Z+ = Viр+, V'(Z) = {dVi/dZj}ij=li2f3. В теореме 1.3.2 формулируются условия при которых решение уравнения (12) можно получить с помощью итерационного процесса Ньютона.

Для дискретизации уравнения (12) по методу Галёркина рассматривается гильбертово пространство Н1(П), состоящее из элементов пространства И^П), для которых выполняется хотя бы одно из условий

j ipdS = 0, или J ipdl = О,

Го Lo

где Го - часть границы Г, a Lo - незамкнутая кривая в области П, со скалярным произведением

(<Р,Ф)& = (Vv?,W0„

Пусть Н^ - m-мерное линейное подпространство Н1 с базисом ..., фт и нормой из Н1. Тогда для аппроксимации Галёркина ¡рт = получим уравнение

Pi&PiVm) - vj) = 0, (15)

где через Рцi обозначен оператор проектирования Н1 на Н^.

Уравнение (15) эквивалентно т-мерной объединённой системе

J kV^pm■Vфidx + J nxVфi J -^—((¡(с^х^ЛЗ^За: = п г г г

г = 1,2,..., ш,

где У(С!)т = {Р(а1)тк — 1)У<рт, а проекция Р(в1)т^'Рт находится из специального уравнения.

Учитывая свойства операторов проектирования, доказываются теоремы 1.3.3 и 1.3.4. Теорема 1.3.3

При выполнении условий теоремы 1.3.1 операторы Рц\Т : Н1 —► Я1 и РТ : Н^ —> Н^ являются сильно монотонными и липшиц-непрерывными. Теорема 1.3.4

При выполнении условий теоремы 1.3.2 производные Фреше РщТ' : Я1 —» Я1 и Р^, Т' :

Н^ —► Н^ являются положительно определёнными и липшиц-непрерывными.

Из теорем 1.3.3 и 1.3.4 следует не только существование и сходимость приближённых решений к точному, но также и сходимость метода последовательных приближений и метода Ньютона.

В § 1.4 рассматривается класс комбинированных систем уравнений для нелинейной задачи (9), полученный с помощью введения вспомогательного вектора. Для простоты изложения предполагается, что граница Г является более гладкой, чем предполагалось раньше, и в каждой точке имеет непрерывный вектор нормали. Для векторов V и 2 имеем

V х 2 = Р, V • V = /, V = кг, х 6 П, (16)

где функция к задана либо как к = к(х,У), либо как к — к(х,2). Используя формулу Грина в области Г2Х для векторов V, 2, получаем

д(х) = -VI ^¿Яу + V х у Зу + V, хеПи (17)

г г

где п - вектор внешней к О нормали. Здесь

П!

Учитывая скачки интегралов в (17) при х —> Г и условия связи нормальной компоненты вектора V и тангенциальной компоненты вектора 2 при переходе через границу Г, приходим к граничному соотношению

= + 1^АЗу + 1[п{У-п)+пх{2хп)), 1бГ. (18)

г г

Здесь

ш = д(х) - \(пМ + Ар х п) + VI ^-¿3у + V х I

г г

аД/ = /-/ь Др = р-р1.

Таким образом, в общем виде задачу (9) можно сформулировать как задачу с уравнениями (16) и граничным условием (18).

С целью получения в (18) только одного интеграла по Г вводится вспомогательный вектор-функция Р, удовлетворяющий в П следующим уравнениям:

V х Р = О, V-P = 0, хеП. (19)

Тогда из формулы Грина имеем

о(х)Р(х) = V х J + V J ^dSy, (20)

г г

где а(х) = 4ж, если х 6 П, и а(х) = 2ж, если х € Г. С учётом (18) и (20), задавая для вектора Р условия на Г в виде Р х п = Z х п или Р • п = V ■ п, получим следующие краевые условия к уравнениям (16) и (19):

Р х п — Z х п, (Р + V ndS) х n = gi(x) x n, (21)

J 4тг r г

P x n = Z x n, + V [ ~ f)' "rfg) • n = й(х) • n, (22)

2 У 47ГГ

г

а также _ _

P -n = V - n, (^4-^ + Vx [(Z~pXndS)xn = g1(x)xn, (23)

2 j 47ГГ

г

P-n = V-H, (P + Vx [ (Z~pXndS)-fi = g1{x)-n, (24)

J 47ГГ г

точно учитывающие поведение решения на бесконечности. Каждая пара краевых условий (21)-(24) содержит только один граничный интегральный оператор и с уравнениями (16), (19) образует четыре различные комбинированные формулировки.

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

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

Как известно, общие подходы к решению нелинейных эллиптических задач были заложены Л.В.Канторовичем, М.М.Вайнбергом, F.E.Browder, G.J.Minty, М.К.Гавуриным. Отметим работы Е.П.Жидкова, И.В.Пузынина, а также М.М.Карчевского, А.Д.Ляшко и других представителей этих школ.

Можно выделить два направления в разработке и применении проекционно-сеточных методов с конечными элементами, которые особенно интенсивно развиваются в последнее время. К первому направлению относятся h—,p— и hp— схемы (р > 2) метода Галёркина с разрывными базисными функциями 6 7, а ко второму - h—,p— и hp— версии схем с векторными конечными элементами ( Whitny's, Nedelec's elements).

"Г.И.Марчук. Методы вычислительной математики. М., Наука, 1977, с. 88-91; 3-е издание, 1989, с.116-122.

'B.Cockburn. Discontinuous Galerkin methods. Z. Angew. Math. Mech. v. 83, N 11, 2003, p.731-754.

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

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

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

2(П) = {и£(И^21(П))3: V • и = О, V х и = 0};

й(П) = {и 6 (12(П))3 : V ■ и = 0, V хие (Ь2(П))3};

У(П) = {и е (Ь2(П))3 : V • и е ¿2(П), V х и = 0}.

Функции из пространства Z(Q) могут быть выражены в виде градиентов гармонических функций. Поэтому в п. 2.1.1 приводятся алгоритмы генерации гармонических конечно-элементных базисных функций.

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

= (25)

¡=1

где и!>, г = 1, ...,т, - значения приближённого решения в узлах .....т на границе ди,

а N¿(1), г = 1, ...,т, - базисные функции.

Для нахождения гармонической базисной функции М(х) такой, что V • УЛ^ = О, согласно первому алгоритму, используется метод коллокации. Представим её в виде

ш ...

х) — Л Чд{])(х)> гДе аУ * неизвестные коэффициенты, д(]) -индексная функция, ¿=1

] = 1,..., т, а в качестве /9(^') выбираются функции

Рп+1л+1(^) = <1пк(г/г0)п со8(к<р)Р*{созв), дп+1М1(х) = (1пк(г/г0)п вт(ку)Р*{созв) (26)

при х = (г, в, уз), (¿„^ = (2п + 1)(п — к)\/(п + к)\, которые вычисляются по реккурентным формулам и входят в общее представление решения задачи Дирихле для оператора Лапласа внутри сферы радиуса г0. Неизвестные коэффициенты находятся в результате решения системы

771

Индексная функция д(]) подбирается так, чтобы система была разрешима. Точность аппроксимаций с помощью N,(1) определяется степенями гармонических многочленов, которые входят в функции /¡,0), ] = 1,..., ш, и для которых аппроксимационная формула (25) точна.

т

По второму алгоритму = ЛЩ'¡¡(х), где выбираются последовательно из

функций (26), а неизвестные коэффициенты находятся в результате решения системы

Е ь1? I = I цью, / = 1,2.....т,

Яш а,.,

ди>

где - лагранжевая граничная базисная функция. Точность приближений функциями ЛЦх) зависит от точности аппроксимаций функциями Li, г = 1 ,...,т.

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

В п. 2.1.2 представлены два алгоритма построения конечно-элементных векторных узловых базисных функции для аппроксимаций высокого порядка из пространств Z(fl)t У(П), С/(Г2). В обоих случаях используется тот факт, что как и сами гармонические функции, их градиенты обладают хорошими аппроксимационными свойствами.

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

Щк)

3 т 3 т

= Е ?*(Е «^ад) = Е Е «у^Л*). (27)

к= 1 ]=1 к=1 ]=1

где и^ - значения приближённого решения в узлах {х3}3=1,...,т на границе ди/.

Согласно первому алгоритму, базисные функции ищуются в виде

где а]*'1' - неизвестные коэффициенты, /9у) - гармонические

функции из (26), д(у) - индексная функция. Неизвестные коэффициенты находятся в результате решения системы

3т г, г Зт д е Зт д л

^<*> =Е^м - о, ЕаГ'&ы'<28>

]=1

цеди, кфкх, кфк2, к\фк2, 1 < к,к\,к2 < 3, / = 1,2,...,т.

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

_ /.V Зт . ..

Во втором алгоритме базисные функции ищуются в виде IV; (г) = ^ Ц V/5у)(х),

где по правилу построения среднеквадратичного приближения коэффициенты б^'1' находятся из системы

3 т

от р *

ви ды

9-^с13, / = 1,2.....Зт. (29)

охк

Здесь д{]) = j + 1, j = 1,2, ...,3т, если Л = 1. Точность аппроксимаций в этом случае зависит от точности приближений лагранжевыми граничными базисными функциями (1 < г < т.).

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

Далее приводятся алгоритмы построения конечно-элементного базиса из гильбертовых пространств С/(П) и Сначала рассматривается случай пространства У(П). Пусть

приближённое решение й ищется в виде (27), где и^ - значения приближённого решения

в узлах {х.,}^.....т внутри и на границе и.

Согласно первому алгоритму, базисные функции представим в виде

¿ а^'''УЛ3ц)(х), где а^'1' - неизвестные коэффициенты, д{]) - индексная функция, Лг(;), ¿=1

] = 1,..., т, - функции из множества

{1; х:; х2; 13; Х1Х2; 12X31X1X2X3; х?, х2, х3,...}.

(30)

Неизвестные коэффициенты находятся в результате решения системы, аналогичной (28).

Зт

Во втором алгоритме базисные функции определяются в виде И7,1*

¿=1

где по правилу построения среднеквадратичного приближения коэффициенты

находятся из системы, аналогичной (29). В обоих алгоритмах точность аппроксимаций

базисными функциями зависит от степени многочленов (30), градиенты которых

участвуют в определении коэффициентов а^''' и

— "* (к) Базисные функции из пространства (/(П) формируются из базисных функций И^ е

К(П), г = 1, ...,т, к = 1,2,3.

В п. 2.1.4 дано более подробное описание конечных элементов с гармоническими

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

элементов, как вид ячейки ш, число узлов т и наибольший порядок (1 гармонических

многочленов, которые точно приближаются соответствующими базисными функциями.

Для стандартного симплекса 5 = {х; > 0, г = 1,2,3,Х1 + х2 + х3 < 1} и пятигранной

призмы Т = {х1, х2 > 0, хг + х2 < 1; -1 < х3 < 1} со стандартным расположением узлов

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

таблицы 2, в случае стандартного куба [—1,1]3 первый алгоритм демонстрирует более

высокую точность.

Таблица 1

5 Т

т 10 20 34 52 18 38

й 2 3 4 5 2 4

Таблица 2

алгоритм 1 алгоритм 2

т 26 56 98 26 56 98

<1 3 5 7 2 3 4

Приведём пример сравнения сходимости приближённых решений в трёхмерной задаче интерполяции магнитного поля, когда в качестве конечных элементов используются последовательности обычных лагранжевых (параллелепипеды с 27, 64, 125 узлами) и гармонических (параллелепипеды с 26,56,98 узлами) элементов. На рис. 1 представлено поведение максимума модуля относительной погрешности й'1' в зависимости от общего числа узлов Л', используемого при интерполяции на адаптивной сетке.

Как видно из рис.1, лучшая точность « 0.1 ■ Ю-4 достигается на 98 узловых гармонических элементах при общем числе точек интерполяции 16358. Это на 37 %

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

Рис. 1. Сходимость интерполяций для параллелепипедов (слева направо, сверху - вниз): о - 27, 64, 125 узловые лагранжевые; • - 26, 56, 98 узловые гармонические (первый метод).

150СЮ 20000 2SOOO

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

V-{a1v) = f, Vxi1=F, i£il; v = 0, хедй.

V -w = /, V x (a2w) = F, x 6 П; w = 0, x e дй.

(31)

(32)

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

Для построения обобщённой формулировки задачи (31) используется вариационный метод. Пусть Й1 = 01(1, и) - заданная функция, / € Ь2(П), Р 6 (Ь2(П)3 и V • Р = 0. Решение задачи ищется в пространстве Ж^Л)3:

Ж,уп)3 = {Йе (12(П))3 : ди/дхк б (¿2(П))3Д = 1,2,3;гГ|вп = 0}.

Как известно, такое пространство позволяет использовать две эквивалентные нормы:

11ж}10(П)3 - llv ' Й11(ЫП))3

+ IIV х ж

11"11^0(п)з = ^2{du/dxk,du/dxk)lL,m з; № к=1

Очевидно, что решение задачи (31) обеспечивает минимум функционала

Fi{v) = (1/2) J [(V ■ (div) - ff + (V x v- Ff}dn.

n

Предположим, что для коэффициента уравнения выполняются следующие условия:

cillV • < ||V • (5iO)||L>,n) < cj|| V • tf||Mn), (33)

|| V • (а^х,?!)«,) - V- (ax(a:, v2)v2)||ia(n> < с3||й - й2||^о(п)з, (34)

где Cj > 0, Ci = const, i = 1,2,3, a v, vi,v2£ W^ffi)3.

18

Справедлива теорема. Теорема 2.2.1

Пусть v+ = tv2 + (1 - t)vi, t е (0,1) и Vi,v2 € W2Если для коэффициента сн выполняются неравенства (33), (34) и справедливо представление

al(x,v+)v+ = tai{x,v2)v2 + (1 - ^(х.й)?!,

то при f 6 L2(U), F € (L2(Cl))3 функционал F\, определённый на И^^П)3, является непрерывным, коэрцитивным и строго выпуклым. При выполнении (33), а также неравенства

|| V ■ (u'(ai,v, и)й)||/,2(п) < йЦиЦи'^щЯ,

где u'{ai,v,u) = {5(oi(i,ff)^)/3iij}jj=i,2,3, v,u£ W^oifO3, c4 = const, доказывается лемма 2.2.1 о том, что градиент функционала Fi, определяемый равенством

(ЛДЙ^одз = J[(V • (ai(x, v)v) - /)V • (u'(auv,u)u) + (V x v- F) ■ V x ujdfi, n

задаёт отображение из Н^0(П)3 в (И^0(П)3)*.

При выполнении условий теоремы 2.2.1 и леммы 2.2.1 из известных теорем о связи метода Рица и метода Галёркина следует, что обобщённая формулировка для задачи (31) в виде уравнения

(AiiT,u) = 0, Vue Wyn)3,

имеет единственное решение v € W^ffi)3. Кроме того, при каждом п это уравнение для Vu„ е (W2i0)„(fi)3 имеет единственное галёркинское приближение йп £ (Н/21о)„(П)3.

Если в задаче (31) F = 0, то для её эффективного решения в полученной обобщённой формулировке рекомендуется использовать базис из пространства V(fi), удовлетворяющий краевым условиям.

Далее вариационный метод используется для построения обобщённой формулировки задачи (32). Предполагается, что а2 = a2(x,w) - заданная функция, / € L2(U), F е (L2(fl))3 и V • F = 0. Решение задачи ищется в пространстве W^0(fl)3. Функционал для задачи (32) имеет вид

Fa(fl) з (1/2) J[(V • (w) - ff + (V х (a2w) - F)2]dil. n

Пусть для коэффициента уравнения выполняются следующие условия:

c5||V х iu|'|(ij(n))3 < ||V х (a2tu)||(Lj(n))3 < c6||V x t5||(bj(n))3, (35)

|| V x (a2(x,uii)wi) - V x (a2(x,iu2)uj2)||(i2(i2))3 < с7\\шг ~ й^Цущз, (36)

где й > 0, с{ = const, г = 5,6,7, a w,wi,w2 е WJi0(fi)3. Справедлива теорема. Теорема 2.2.2

Пусть w+ = tw2 + (1 - t)wi, t 6 (0,1) и щ,ги2 е Ж^П)3. Если для коэффициента а2 выполняются неравенства (35), (36) и справедливо представление

a2(x,w+)w+ = ta2{x,w2)b}2 + (1 - t)a2(x, щ)^,

то при J 6 L2{Q.), F £ (L2(Q))3 функционал F2, определённый на Wj^fi)3, является непрерывным, коэрцитивным и строго выпуклым.

Обобщённая формулировка для задачи (32) получается с помощью дифференцирования функционала F2. При выполнении (35) и неравенства

IIV х (v'(a2,w, u)u)||i,2(n) < свЦйЦ^^з,

где v'(a2,w,u) = {d{a2(x,w)wi)/duj}ij=w, w,u € W2,o(n)3, с8 = const, доказывается лемма 2.2.2 о том, что градиент функционала F2, определяемый равенством

(A2w,u)(L2m3 = J|(V ■ w - /)V • й + (V x (a2{x,w)w) - F) • V x (i/(52,ii;,u)u)]dfi, n

задаёт отображение из H^offi)3 в (И^0(П)3)*.

При выполнении условий теоремы 2.2.2 и леммы 2.2.2 обобщённая формулировка для задачи (32) в виде уравнения

(A2w, й) = о, Weiy2yn)3,

имеет единственное решение iu е W^f})3. Кроме того, при каждом п это уравнение для Vu„ £ (H/2io)n(ii)3 имеет единственное галёркинское приближение wn 6 (W2o)„(fi)3.

Если в задаче (32) / = 0, то для её эффективного решения в полученной обобщённой формулировке_рекомендуется использовать удовлетворяющий краевым условиям базис из пространства U(П).

В § 2.3 представлены проекционно-сеточные схемы для решения нелинейных задач вида (31),(32) относительно вектор-функций, использующие предложенные в п.2.1.3 функции формы из пространства К(П). Основой для построения этих схем является подход Бубнова-Галёркина.

Здесь строится проекционно-сеточная схема относительно вектор-функции для решения нелинейной краевой задачи с внутренней границей:

V ■ V = О, V х Z = 0, V = pZ, х е П;

V-V = 0, Vx2 = F' (z), V = Z, ieil'; (37)

[n • V\ = 0, [n x Z\ = 0, х£дП; V = Z = 0, x £ Г0.

Предполагается, что граница дП является общей для областей f! и П', причём Г)' окружает П так, что Г0 является удалённой от Г2 частью границы 8Q'. Важно отметить, что область Q может быть составлена из нескольких несвязных областей. Функция р задана либо как р = p(x,V), либо как р = p{x,Z). Через [п • V] и [n х Z] обозначена разность пределов соответствующих функций при стремлении точки х на границу 8Q изнутри и извне области П. Относительно нелинейного коэффициента задачи предполагается, что выполняются обычные неравенства.

Дискретизованная обобщённая формулировка задачи (37) обеспечивает выполнение условий на внутренней границе для нормальных _и тангенциальных компонент искомой вектор-функции. С базисом из пространства V(U) при формировании матрицы не требуется вычисления слагаемых, содержащих ротор.

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

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

Третья глава диссертации посвящена методам и алгоритмам, существенно повышающим эффективность решения как прямых, так и обратных нелинейных задач магнитостатики. Среди публикаций, посвящённых построению проекционно-сеточных схем для задач магнитостатики, отметим работы Н.И.Дойникова, С.Е.Сычевского, Е.А.Ламзина с соавторами, J.C.Coulomb с соавторами, J.Simkin, C.W.Trowbridge с соавторами, M.J.Friedman, J.E.Pasciak, I.D.Mayergoyz с соавторами, Е.П.Жидкова, С.Б.Ворожцова, П.Г.Акишина, Б.Н.Хоромского, Э.А.Айряна с соавторами, А.Г.Дайковского с соавторами, S. Russenschuck с соавторами и других исследователей.

В § 3.1 предлагается новый метод вычисления функции на липшицевой границе трёхмерного тела по заданному градиенту. В качестве примера применения метода приведём следующую нелинейную краевую задачу. Пусть область П' С R3 ограничена липшицевой границей и П' = fílJOi, где области Q и Q¡ разделяются внутренней границей Г С П', которая также является липшицевой. Требуется найти функции и, щ, удовлетворяющие уравнениям

3 д

У —(аЛх, и, Vu)) + olí, u, Vu) = 0, х ей;

, OXi

1=1

+ *6П1; (38)

ij=l J i= 1

з 3 du

и — txi = 0, ^2{ai(xlu,'Vu)— ^2äij(x)-~ — biUi) ■ a¡ = 0, х e Г; ¡=i j=i öxi

3 du 3

u = ü(x), x e <ЭП \ Г; (^ + ^Mi) ■ 4 = üi(i), x6öfii\r,

i ÖXj .

1,3 = 1 J 1=1

где Vu - градиент функции и, а aj,á¿ - компоненты векторов внешних нормалей к соответствующим областям. При этом равенства на границе понимаются в смысле следов соответствующих функций. Предполагается, что все коэффициенты в уравнениях вещественные и вместе с заданными функциями /,й,щ удовлетворяют условиям, необходимым для однозначной разрешимости приведённой задачи. Если в области fii известно частное решение - достаточно гладкая вектор-функция G, компоненты которой Gi (г = 1,2,3) удовлетворяют уравнению

3 д 3

£ дГЫ*) ' Gj) + £ bí ■ Gi = /, X 6 пь

i,j=1 ' t=l

то можно ввести функцию ф, такую, что Уф = G в ÍV Эта функция будет являться частным решением уравнения в области fii. Тогда задачу (38) можно сформулировать

относительно функций и в области Q и йх = щ — ф в области П]. Причём в Г2Х получается однородное уравнение для функции щ. Поэтому если вектор G вычисляется более точно и быстро, чем функция / в области fii, и при этом упрощается процесс построения трёхмерной сетки, то переформулировка задачи для неизвестных и,щ является оправданной. Однако точность решения новой задачи зависит от точности вычисления функции ф, которая входит в условия на границе Г. Для нахождения функции ф вводится пространство функций

IV2(Г) = {и : и 6 ¿2(Г), VTu 6 (i2(r))2},

где в локальной правой декартовой системе координат с ортами s, t, п, направленными вдоль касательной, бинормали и внешней нормали,

VTu = sdu/ds + tdu/dt.

В Ж2' (Г) выделяется класс функций И^'оСГ), для которых выполняется хотя бы одно из условий

J udS = 0 или J udl = О,

Го !о

здесь Г0 - часть границы Г, а io - некоторая незамкнутая кривая на Г. Для функций u,v 6 Ж20(Г) можно ввести скалярное произведение и соответствующую норму

(u,v)wio( Г) = (VTu, VTu)(i,2(r))2, |Mliy2io(r) = ((и.иЦуг))172-

Введённая норма эквивалентна норме пространства Ж2 (Г). Функция ф находится из обобщённого уравнения

j Угф ■ VTvdS + J p(G, |VT0|)(VT0 - GT) • VTvdS = J Gr ■ VTvdS, v e Ж2уг),

Г1 Г2 Г,

где Гь Г2 - непересекающиеся составные части границы Г, определяемые в зависимости от поведения G е К = {Уф \ ф е Ь2(Г),Чф 6 (£2(Г))3)}. Здесь

p(G, |VT0|) = 2 |Gn|

|G-VT0|

Для исследования этого уравнения вводятся следующие обозначения. Зададим операторы В], В2 с помощью равенств

(Вхи^Цуг,) = J ■ ЧгчЛБ, 6 И^>0(Г1);

Г1

(В2и2,^)и,.о(г2) = I р{б, |Ути2|)(Ути2 - <3Т) ■ УтУ^Б, и2,ь2 £ ЖзУГг). г2

Очевидно, что В, : Ж^Г^) -» Ж^) и В2 : И^о(Га) -» Ж2уг2). Свойства оператора В2 сформулированы в следующей лемме. Лемма 3.1.1 Оператор В2 : Ж20(Г2) —» Ж210(Г2) является:

1) сильно монотонным, то есть выполняется неравенство

(В2и - BiV,U- 1))wl0(r2) > II" - 'ullw2yr2)> N2,0^2);

2) липшиц-непрерывным, то есть выполняется неравенство

||В2и - B2u||wio(rj) < <т7||и - w||w2yra). и, v 6 И^10(Г2), где сг7 = 2 + max(l/|G„|), бе К.

хеГг

На основе леммы 3.1.1 доказывается, что оператор В : Я —» Я, где Я = (П ni^i^i) х ^2(Г2)}, задаваемый равенством

(.BU,V)H = (В^.и^уг,) + (В2И2,^)и/210(Г2)

для u = (ui,u2), v = (vi,v2) е Н, является сильно монотонным и липшиц-непрерывным. Из общих теорем об уравнениях с такими операторами сразу же получаем однозначную разрешимость обобщённого уравнения задачи и сходимость галёркинских приближений к точному решению. Кроме того, известен итерационный процесс, с помощью которого можно получить это решение.

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

В этом же параграфе приводится алгоритм вычисления функции на липшицевых границах несвязных областей по заданному градиенту для задач магнитостатики. Скалярный потенциал проводников с током <ps предлагается определять на границах Yk,k — 1,2,..., т, несвязных областей из уравнений

Vy>s = Hs{х), хеГк, к = 1,2, ...,т, (39)

с условием <ps{yk) = vfi Ук £ Г*, tps(yo) - известная величина, где уо может не принадлежать границам = 1,2, ...,tti. В общем случае точки ук,к = 1,2, ...,т,

связаны с точкой у0 ломаными Ьк,к = 1,2,...,т, а потенциал ¡ps на этих ломаных

находится с помощью адаптивного алгоритма из обобщённых уравнений

.....

U Lk

где v - достаточно гладкая функция, Хк - функция, обеспечивающая выполнение условия в точке я/о> вектор t(i) задаёт направление ломаной Обобщённые линейные уравнения, соответствующие (39), будут иметь вид

J VT/ • VTvdS = J(Hs + VTXk) • VTvdS, к = 1,2,..., m, (40)

rt rt

8И.Н.Бронштейн, К.С.Семендяев.Справочник по математике,10-е издание,М.: Наука, 1965, с.538-539.

где функции v, Хк имеют тот же смысл, что иг;, \к в предыдущем уравнении соответственно.

Для решения уравнений (40) адаптивным методом конечных элементов рассматриваются четыре типа четырехугольных элементов. При построении конечно-элементной системы используются лагранжевы элементы с 9-ю узлами (тип 2) и серендиповы элементы с 17-ю узлами (тип 3). Для вспомогательных вычислений применяются лагранжевы элементы с 6-ю (тип 1) и 15-ю (тип 4) узлами. Доказана лемма 3.1.3 об оценках функции интерполянта через старшие производные известной вектор-функции. Эффективность адаптивного алгоритма демонстрируется на модельном примере и на вычислении скалярного потенциала обмотки для большого соленоидального магнита L3 [1].

§ 3.2 диссертации посвящен алгоритмам с контролем точности для решения задач магнитостатики. Разработки алгоритмов контроля точности для проекционно-сеточных методов содержатся в многочисленных работах. Особенно актуальным этот вопрос является для численного решения трёхмерных эллиптических задач с сильно меняющимися решениями, в расчётной области которых содержатся подобласти с неодинаковыми дифференциальными уравнениями. При этом сами подобласти могут иметь сложную конфигурацию. Здесь можно отметить книги С.Г.Михлина, Х.Л.Смолицкого(1965), Н.С.Бахвалова(1975), Н.Н.Калиткина(1978), Г.И.Марчука,В.В.Шайдурова(1979), В.В.Шайдурова(1989), B.Szabo, I.Babu8ka(1991), Ch.Schwab(1998), O.C.Zienkiewicz, R.T.Taylor(2000), P.Solin(2006), в которых рассматриваются вопросы, возникающие при построении такого типа алгоритмов и есть ссылки на оригинальные публикации. Различные алгоритмы контроля точности используются в некоторых коммерческих системах автоматизированного проектирования при решении линейных и нелинейных краевых задач с эллиптическими дифференциальными операторами второго порядка. Все апостериорные методы анализа ошибки проекционно-сеточных методов можно разделить на три основные группы: 1) методы, основанные на решении задачи в дополняющих друг друга формулировках, когда точное решение принимает промежуточное значение между приближёнными решениями одной и другой формулировки; 2) методы с многократным решением задачи; 3) методы с локальным оцениванием ошибки. Наиболее экономичной по вычислительным затратам считается третья группа методов, которая и развивается в последнее время особенно интенсивно. В этой группе выделяются методы с целенаправленной /ip-адаптивностью (L.Demkowicz с соавторами), а также метод суперсходящегося покусочного восстановления (superconvergence patch recovery method, O.C.Zienkiewicz, J.Z.Zhu). Об актуальности вопросов, рассматриваемых в этой главе, свидетельствует практически постоянное представление докладов по контролю точности расчётов в задачах магнитостатики на регулярных конференциях Compumag, посвящённых моделированию магнитных полей и смежным темам. В отличие от известных методов для этого класса задач, в третьей главе предлагается комплексный подход, объединяющий на основе единых характеристик как тестирование функций магнитного поля, основанных на данных измерений для спектрометров физических экспериментов, так и тестирование магнитного поля, которое получается в результате расчётов.

В параграфе приводятся формулы, характеризующие численную погрешность в зависимости от невязки уравнений для операторов дивергенции и ротора. Пусть область U С R3 разбита стандартным образом на некоторые конечные элементы W{, г — 1,2,..., N, и на каждом г-ом элементе определены функции формы Nitk(x), так что N^ix^) = ¿(x^xW) где узлы, х^ух^ е Тй<. Через wk обозначим носитель Ахгой базисной

функции тк = и виррИ^к, а через (Со°(й/к))3 - множество финитных бесконечно-

дифференцируемых вектор-функций с носителями на йк. Справедлива следующая лемма. Лемма 3.2.1

Пусть для V 6 (Ср0^))3 выполняются неравенства

IV • ф) - V • %)| < схйо, |У х ф) - V х %)| < с2с!о, где х, у 6 ¿о = |х — г|, £1, с2 - положительные постоянные, не зависящие от

х и у. Тогда при х € шк выполняется равенство

ф) = (V ■ ф))/0(х) + (V х ф)) х Га{х) + До?,

здесь

То{х) = ^ I ^¿у, |/0(х)| < (¡о, < с34,

¿к 3 _

при этом с3 = тах{с1,с2}, г = £ - УР), г = \х - у|.

Р=1

Справедлива также лемма о представлении вектор-функции V е (С°°(и^))3 на элементе

и>к-

Лемма 3.2.2

Пусть для V 6 (С°°(ш(ь))3 выполняются неравенства

IV • ф) - V ■ %)| < с4<*, IV х ф) - V х ф)| < с5й, где х, у £ ик, й — тах |х — г\, £4, сд - положительные постоянные, не зависящие от х и у. Тогда при х 6 и>к выполняется равенство

ф) = (V ■ ф))Т(х) + (V х ф)) х Т(х) + Ру + №,

здесь

_1 г - 3

/(1) = / 1/(1)1 - - -

и,к У Р=!

при этом 56 = тах{с4, £5}.

Леммы 3.2.1 и 3.2.2 можно использовать для оценки погрешности приближённого решения краевых задач с операторами дивергенции и ротора.

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

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

т}(/("в,ш) = (| у V ■ /(,)в<н +11V х /<'>В<Н)/(МВ0),

IV Ъ)

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

быть равен оператору 1е \ задаваемому формулой ii^B = ¿ДЛ^(х), где индекс I

¿=1

- степень интерполяционного многочлена, mi - число используемых при интерполяции узлов сетки, Д - значение вектора В в i-м узле, N}l\x) - г-я базисная функция интерполяционного многочлена.

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

т;(/е(,)Д w) < С ■ hM ■ \DMB{x)\, few,

где С = const, С > 0, h - диаметр элемента w, Dm = ö|m|/dmixi&П2х2дтзх3, \т\ = тх + т2 + т3.

В § 3.3 представлены методы решения задачи формирования однородного магнитного поля, которую можно отнести к классу обратных задач. Общий подход к их решению разработан А.Н.Тихоновым и его школой. Здесь следует отметить работы П.Н.Вабищевича, а для задач магнитостатики - работы Е.А.Ламзина, С.Е.Сычевского с соавторами. В этом разделе сформулирован итерационный процесс уточнения части поверхности ферромагнетика для получения высокооднородного магнитного поля (теорема 3.3.1) и даны примеры решения задач с его помощью.

В § 3.4 приводится алгоритм построения компьютерных моделей анализирующих магнитов спектрометров с контролем точности вычислений. Главными требованиями к таким магнитам обычно являются: величина поворотной силы, объём рабочей области, эксплуатационные условия, материал проводника и тип стали. Процесс построения компьютерной модели дипольного магнита, основанный на этих входных данных, можно описать с помощью следующего алгоритма:

1) оценка ампер-витков, длины магнита и толщины ярма по аналитическим формулам;

2) построение компьютерной модели обмотки и итерационный процесс формирования геометрии ненасыщенного магнитного ярма;

3) анализ точности разработанной компьютерной модели, который включает в себя проверку правильности вычисления поля от обмотки и анализ адекватности расчётной сетки;

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

Четвертая глава диссертации посвящена результатам применения методов, описанных в предыдущей главе, к построению компьютерных моделей и расчётов полей магнитных систем для установки ALICE [1-4] (CERN, г. Женева) и проекта установки для эксперимента с поляризованной мишенью в ИТЭФ [5] (г. Москва).

В § 4.1 приводятся результаты построения компьютерных моделей для большого соленоидального магнита L3, проекта сверхпроводящего дипольного магнита и первого проекта "теплого" дипольного магнита для эксперимента ALICE. В частности, для модели

L3 представлены результаты сравнения расчётов, выполненных по программе MSFE3D 9, с картой поля магнита, а также с результатами расчётов по другим программам.

В § 4.2 представлены результаты построения компьютерных моделей дипольных магнитов серии VULCAN с конусной конфигурацией ферромагнетика для эксперимента ALICE. Здесь приводятся разработанные модели обмоток, укладываемых на боковую поверхность прямого усечённого конуса, результаты использования различных характеристик для оценки точности вычисления поля в апертурах магнитов, а также другие расчёты.

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

К диссертационной работе имеется три приложения. Приложение 1 содержит анализ функций магнитного поля, основанных на данных измерений для эксперимента ЭКСЧАРМ [6-8]. Рассматриваются 4 различные функции магнитного поля, интерполирующие данные с сетки измерений с помощью произведения интерполяционных многочленов Лагранжа первой, второй третьей и четвертой степеней соответственно по каждой пространственной переменной. Для их сравнения используется характеристика Г){№В, w0), 1 < к < 4, введённая в § 3.2.

В Приложении 2 приводятся результаты расчётов магнитного поля для модели соленоида проекта модифицированного бетатрона, разработанного в Лаборатории Ядерных Проблем ОИЯИ.

В Приложении 3 представлены результаты расчётов магнитного поля для модели поворотного магнита М-90-05 ускорительного комплекса У400-У400М, предназначенного для проекта DRIBs [9].

ОСНОВНЫЕ РЕЗУЛЬТАТЫ ДИССЕРТАЦИИ

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

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

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

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

3. Разработаны методы и алгоритмы решения трёхмерных нелинейных задач магнитостатики:

- проекционно-сеточный метод и адаптивный алгоритм вычисления скалярного потенциала, создаваемого проводниками с постоянным током, на границах несвязных областей вне проводников;

- алгоритм построения единых характеристик для контроля точности вычисленных и измеренных магнитных полей ;

- итерационный метод формирования высокооднородного магнитного поля за счёт изменения объема, занимаемого ферромагнетиком;

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

4. На основе созданных алгоритмов и программного обеспечения решён ряд трёхмерных задач с учётом нелинейного влияния ферромагнитных частей магнитных систем. В некоторых задачах использовались характеристики точности вычислений. Разработаны компьютерные модели и проведены расчёты:

- большого соленоидального магнита L3 установки ALICE на LHC (CERN, г.Женева) и проектов ОИЯИ дипольных магнитов для этой установки: сверхпроводящего и "тёплых" вариантов, в том числе с конусной конфигурацией магнитопровода;

- спектрометрического магнита с поляризующими наконечниками для проекта эксперимента с поляризованной мишенью на ускорителе в ИТЭФ (г.Москва);

- соленоида для проекта модифицированного бетатрона ОИЯИ (г.Дубна).

Диссертация основывается на следующих опубликованных работах

1. Е.П.Жидков, А.В.Фёдоров, О.И.Юлдашев. Об одном подходе к решению задачи маг-

нитостатики в комбинированной постановке. Математическое моделирование, т.2, N 9, 1990, с. 10-20.

2. Е.П.Жидков, А.В.Фёдоров, О.И.Юлдашев. Использование метода Галёркина для

решения интегрального уравнения магнитостатики относительно векторного потенциала. Дифференциальные уравнения, т.28, N 10, 1992, с. 1821-1828.

3. Е.П.Жидков, О.И.Юлдашев, М.Б.Юлдашева. Адаптивный алгоритм вычисления

функции на липшицевой границе трёхмерного тела по заданному градиенту и его применение в магнитостатике. Журнал вычислительной математики и математической физики, т.42, N 12, 2002, с.1854-1869.

4. E.P.Zhidkov, O.I.Yuldashev, M.B.Yuldasheva. A projection method for solving linear prob-

lems with the divergence, curl operators and its application in magnetostatics. Вестник Российского университета дружбы народов. Серия Прикладная и компьютерная математика, N 1(1), 2002, с.79-86.

5. Е.П.Жидков, О.И.Юлдашев, М.Б.Юлдашева. Новые проекционные формулировки

относительно векторов поля для решения нелинейных задач магнитостатики. Вестник Российского университета дружбы народов. Серия Прикладная и компьютерная математика, т.2, N 2, 2003, с. 104-115.

6. Е.П.Жидков, В.В.Рыльцов, О.И.Юлдашев, М.Б.Юлдашева. Численное решение зада-

чи формирования однородного магнитного поля за счёт изменения занимаемого ферромагнетиком объёма для некоторых магнитных систем экспериментальной физики. Вестник Российского университета дружбы народов. Серия Физика, N 12, 2004, с.17-25.

7. Е.П.Жидков, О.И.Юлдашев, М.Б.Юлдашева. О контроле точности вычислений при

моделировании пространственных магнитных полей. Вестник Российского университета дружбы народов. Серия Прикладная и компьютерная математика, т.4, N 1, 2005, с.93-101.

8. О.И.Юлдашев, М.Б.Юлдашева. Об одном классе конечных элементов с гармоничес-

кими базисными функциями. Вестник Российского университета дружбы народов. Серия Математика. Информатика. Физика, N 2 (2), 2010, с.44-48.

9. A.Ivanov, S.Ivashkevich, ... I.Meshkov, ... O.Yuldashev. Focusing system of the modified

betatron: design, technology, manufacturing and test. Nuclear Instruments and Methods in Physics Research, section A, v.441, 2000, p.262-266.

10. P.G. Akishin, E.V. Arkhipov, ... A.S.Vodopianov, ... O.I. Yuldashev. Superconducting

dipole magnet for ALICE dimuon arm spectrometer. Краткие сообщения ОИЯИ, N 1, 1997, p.81-94.

11. Э.А.Айрян, Е.П.Жидков, Б.Н.Хоромский, О.И.Юлдашев. Алгоритм учёта условий

на бесконечности в двумерных задачах магнитостатики. ОИЯИ, PI 1-87-49, Дубна, 1987, 15с.

12. Е.П.Жидков, В.В.Журавлев, В.С.Кладницкий, И.М.Матора, A.B. Фёдоров, О.И.Юл-

дашев. О формировании однородного магнитного поля на инжекторном участке ЛИУ-30. ОИЯИ, Р9-88-508, Дубна, 1988, 11с.; http://ccdb4fs.kek.jp/cgi-bin/img_index?8811447.

13. Е.П.Жидков, А.В.Фёдоров, О.И.Юлдашев. Алгоритм численного моделирования

фокусирующей магнитной системы ЛИУ. Тезисы докладов Всесоюзной конференции "Вычислительная физика и математическое моделирование", Волгоград-1989, УДН, М., 1989, с.46-47.

14. Е.П.Жидков, В.В.Журавлев, В.С.Кладницкий, И.М.Матора, A.B. Фёдоров, О.И.Юл-

дашев. О формировании однородного магнитного поля на инжекторном участке ЛИУ-30 с использованием ферромагнитных колец. Труды XI Всесоюзного совещания по ускорителям заряженных частиц. ОИЯИ, Д9-89-52, Дубна, 1989, т.1, с.381-383.

15. E.A.Ayrjan, E.P.Zhidkov, A.V.Fedorov, O.I.Yuldashev. A Galerkin Method for Solving

the Magnetostatic Nonlinear Potential Integral Equations. Proceedings of International Symposium on Numerical Análisis (ISNA'92). Part III, Contributed Papers. Prague, August 31-September 4, 1992, p.20-32.

16. E.A.Ayrjan, A.V.Fedorov, O.I.Yuldashev, E.P.Zhidkov. One approach for coupling FEM

and BEM. Proceedings of International Symposium on Numerical Análisis (ISNA'92). Part III, Contributed Papers. Prague, August 31-September 4, 1992, p.12-19.

17. P.G.Akishin, A.S.Vodopianov, I.V.Puzynin, Yu.A.Shishov, M.B.Yuldasheva, O.I.Yulda-

shev. Computing Models of the L3 Magnet and Dipole Magnet for the ALICE Experiment. CERN-ALICE 96-06, Int. Note/Mag, 4 April, 1996, 21p.

18. P.G. Akishin, E.A. Arkhipov, ... A.S.Vodopianov, ... O.I. Yuldashev. Conceptual design

of the warm dipole magnet for the ALICE forward muon spectrometer. CERN-ALICE 96-26, Int. Note/Mag, 1996, Юр.

19. П.Г.Акишин, А.С.Водопьянов, И.В.Пузынин, Ю.А.Шишов, М.Б.Юлдашева, О.И.Юл-

дашев. Моделирование поля мюонного детектора установки ALICE. XV Совещание по ускорителям заряженных частиц, сборник докладов. ГНЦ РФ ИФВЭ, Протвино, 1996, т.2, с.182-185.

20. A.S. Vodopianov, Yu.A. Shishov, М.В. Yuldasheva, O.I. Yuldashev. Computer models of

dipole magnets of a series "VULCAN" for the ALICE experiment. JINR, Ell-98-385, Dubna, 1998. 18p.; http://cdsweb.cern.ch/search?sysno=000328897CER

21. Р.И.Давлетшин, Е.П.Жидков, В.В.Куликов, В.В.Рыльцов, О.И.Юлдашев, М.Б.Юл-

дашева. Компьютерная модель магнита для проекта эксперимента с поляризованной мишенью в ИТЭФ. Расчёт основных вариантов без поляризующих наконечников. ОИЯИ, Р11-98-351, Дубна, 1998, 18с.

22. J.Ritman, O.I.Yuldashev, M.B.Yuldasheva. An algorithm for construction of dipole mag-

nets computer models with quality control and its application for the PANDA Forward Spectrometer. JINR, Ell-2005-49, Dubna, 2005, 19p.

23. О.И.Юлдашев, М.Б.Юлдашева. Новый проекционно-сеточный подход для моделиро-

вания пространственных нелинейных магнитных полей сложных магнитных систем экспериментальной физики. Международная конференция "Тихонов и современная математика", Тезисы докладов секции Математическое моделирование, с.189-190. М., МГУ, 2006.

24. О.И.Юлдашев, М.Б.Юлдашева. Гармонические базисные функции для конечных

элементов высокого порядка аппроксимации. JINR LIT Scientific Report 2006-2007, JINR, Dubna, 2007, c.317-320.

25. О.И.Юлдашев, М.Б.Юлдашева. О конечно-элементном подходе относительно век-

торов поля для расчётов сложных магнитных систем экспериментальной физики. JINR LIT Scientific Report 2006-2007, JINR, Dubna, 2007, c.234-238.

26. O.I.Yuldashev, M.B.Yuldasheva. 3D finite elements with harmonic basis functions for

approximations of high order. JINR, El 1-2008-104, Dubna, 2008, 24p.

27. О.И.Юлдашев, М.Б.Юлдашева. Новое правило построения квадратичных функ-

ционалов для решения методом конечных элементов относительно вектор-функций краевых задач с системой из двух уравнений первого порядка в дивергентной и вихревой формах. JINR LIT Scientific Report 2008-2009, JINR, Dubna, 2009, c.113-116.

28. О.И.Юлдашев, М.Б.Юлдашева. Конечно-элементные векторные узловые базисные

функции из специальных гильбертовых пространств. JINR LIT Scientific Report 2008-2009, JINR, Dubna, 2009, c. 105-108.

Цитируемая литература

1. ALICE collaboration (N.Ahmad, ... O.I.Yuldashev et al.) ALICE Technical Proposal,

CERN/LHCC 95-71, 1995, 237p.

2. ALICE collaboration (S.Beole, ... O.I.Iouldachev et al.) The Forward Muon Spectrometer

of ALICE (Addendum to the ALICE technical proposal), CERN/LHCC 96-32, 1996, 77p.

3. ALICE collaboration. (S.Beole, ... O.I.Iouldachev et al.) ALICE technical design report:

detector for high momentum PID. CERN-LHCC-98-19, 1998, 198p.

4. ALICE collaboration. (G.Dellacasa, ... O.I.Iouldachev et al.) ALICE technical design

report: photon multiplicity detector (PMD). CERN-LHCC-99-32, 1999, 134p.

5. И.Г.Алексеев, Н.А.Бажанов, ... В.В.Рыльцов и др. Исследование поляризационных

параметров в бинарных реакциях рождения странных частиц пр —> КА(Е) на ускорителе ИТЭФ (предложение эксперимента). ИТЭФ 41-97, Москва, 1997, 33с.

6. ЭКСЧАРМ коллаборация ( А.Н.Алеев, ... О.И.Юлдашев и др.). Исследование инклю-

зивного образования ip - мезонов нейтронами на серпуховском ускорителе. ОИЯИ, Р1-96-437, Дубна, 1996.

7. ЭКСЧАРМ коллаборация ( А.Н.Алеев, ... О.И.Юлдашев и др.). Наблюдение очаро-

ванного бариона в эксперименте ЭКСЧАРМ. Краткие сообщения ОИЯИ, N 3 [77], 1996, с.31-46.

8. ЭКСЧАРМ коллаборация ( А.Н.Алеев, ... О.И.Юлдашев и др.). Спектрометр

ЭКСЧАРМ. ПТЭ, N 4, 1999, с.52-64.

9. G.G.Gulbekian, Yu.Ts.Oganessian. "DRIBs": The Dubna Project for Radiactive Ion Beams.

Nuclear Shell - 50 Years: Intern. Conf. on Nuclear Physics, Dubna, Apr. 1999. Proc.-Singapore etc.:World Sci., 2000.

Получено И августа 2010 г.

Отпечатано методом прямого репродуцирования с оригинала, предоставленного автором.

Подписано в печать 13.08.2010. Формат 60 х 90/16. Бумага офсетная. Печать офсетная. Усл. печ. л. 2,12. Уч.-изд. л. 3,57. Тираж 100 экз. Заказ № 57064.

Издательский отдел Объединенного института ядерных исследований 141980, г. Дубна, Московская обл., ул. Жолио-Кюри, 6. E-mail: publish@jinr.ru www.jinr.ru/publish/

Введение

Глава 1. Проекционно-сеточные методы решения краевых задач для векторов с интегральными граничными условиями, точно учитывающими убывание решения на бесконечности

1.1 Теорема о компонентах непрерывно-дифференцируемого вектора и её следствия.

1.2 Интегральное уравнение для векторного потенциала.

1.2.1. Исследование уравнения (1.2.1) в пространстве Н(го1;,П)

1.2.2. Решение уравнения (1.2.1) методом Галёркина.

1.3 Интегральное уравнение для скалярного потенциала.

1.3.1. Свойства интегрального оператора и производной Фреше для скалярного потенциала.

1.3.2. Метод Галёркина для уравнения (1.3.1).

1.4 Комбинированные формулировки для задачи (1.1.9).

1.4.1. Комбинированные формулировки и введение дополнительного вектора.

1.4.2. Проекционная формулировка задачи (1.4.10),(1.4.13)

1.4.3. Пример расчёта плоскопараллельных магнитных полей с помощью комбинированной формулировки.

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

2.1 Векторные узловые конечно-элементные базисные функции из специальных гильбертовых пространств.

2.1.1. Общие свойства гармонических конечно-элементных базисных функций.

2.1.2. Гармонические конечные элементы в R3.

2.1.3. Конечно-элементные базисные функции из пространств Z{ÇÏ), Û{Cl).

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

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

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

2.2.1. Система с нелинейным уравнением в дивергентной форме

2.2.2. Система с нелинейным уравнением в вихревой форме

2.3 Проекционно-сеточные схемы относительно вектор-функций, использующие базис из пространств V(f2), U(Q).

2.3.1. Общие свойства проекционно-сеточных схем, использующих базис из пространства V(£î).

2.3.2. Проекционно-сеточная схема для решения нелинейной краевой задаче с внутренней границей.

2.3.3. Решение относительно векторного потенциала нелинейной краевой задачи с внутренней границей.

Глава 3. Алгоритмы с контролем точности вычислений для решения задач магнитостатики.

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

3.1.1. Проекционно-сеточный метод.

3.1.2. Адаптивный алгоритм для решения системы (3.1.12)

3.1.3. Применение проекционного метода и адаптивного алгоритма в магнитостатике.

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

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

3.2.2. Характеристика функции магнитного поля, основанной на данных измерений.

3.3 Методы решения задачи формирования однородного магнитного поля.

3.3.1. Итерационный процесс уточнения поверхности ферромагнетика для получения высокооднородного магнитного поля. •

3.3.2. Повышение однородности магнитного поля на примере линейного индукционного ускорителя.

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

3.4.1. Аналитическая оценка ампер-витков, длины магнита и толщины ярма.

3.4.2. Построение компьютерной модели обмотки и формирование ярма магнита.

3.4.3. Контроль точности компьютерной модели.

3.4.4. Вывод основных характеристик магнита.

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

4.1 Построение моделей для магнита L3, проекта сверхпроводящего дипольного магнита и первого проекта теплого дипольного магнита для эксперимента ALICE.

4.2 Компьютерные модели проектов дипольных магнитов серии VULCAN для эксперимента ALICE. *

4.2.1. Построение компьютерных моделей обмоток.

4.2.2. Построение компьютерных моделей магнитов с локальным контролем точности.

4.2.3. Результаты расчётов.

4.3 Компьютерная модель магнитной системы для проекта эксперимента с поляризованной мишенью в ИТЭФ

4.3.1. Моделирование поля проекта магнита без поляризующих наконечников.

4.3.2. Расчёты сил и оптимизация поверхности поляризующих наконечников.

Введение диссертационной работы состоит из пяти разделов. Первый посвящён актуальности рассматриваемых проблем. Во втором разделе приводятся основные утверждения общей теории проекционно-сеточных методов для эллиптических задач с монотонными операторами, которые используются в диссертации. В третьем даётся краткий обзор методов для решения эллиптических задач с контролем точности. В четвертом разделе представлен краткий обзор современных подходов к расчётам стационарных магнитных полей для установок экспериментальной физики. Наконец, в пятом приводится краткое содержание диссертации по главам, проведённая апробация, публикации, на основе которых написана диссертация и принятая в работе нумерация. Во введении и далее используются следующие эквивалентные обозначения для дифференциальных операторов векторного анализа: grad а = Va, div{grad а) = V2a - для скалярной функции a; div а = V - о, rot a = Vxa, grad(div а) — rot(rot а) = V2a - для вектор-функции а.

1. Актуальность проблем, рассматриваемых в диссертации.

Проекционно-сеточные методы получили широкое распространение в практических приложениях, благодаря их успешному применению в области теории упругости в начале прошлого века. В настоящее время этот класс методов является наиболее распространённым в инженерных расчётах и с успехом используется в самых различных областях теории и практики. В связи с этим следует упомянуть А.Н.Крылова, Л.В.Канторовича, В.И.Крылова, И.С.Березина, Н.П.Жидкова, Г.И.Марчука, A.A.Самарского, А.В.Гулина, Н.С.Бахвалова - авторов основополагающих книг по численным методам [1]-[7] и других.

Отметим наиболее общие достижения за последние два десятилетия в развитии проекционно-сеточных методов для решения эллиптических краевых задач: это различные алгоритмы решения дискретизованных нелинейных уравнений на последовательности сеток; р— и hp— версии метода Бубнова-Галёркина для непрерывных и разрывных базисных функций, а также для векторных рёберных и граневых (facet) базисных функций; кроме того, это смешанные методы, методы декомпозиции области, методы расчёта с использованием комбинированных систем (т.е. систем с дифференциальными и граничными интегральными уравнениями). Все эти подходы к численному решению задач в различных областях знаний интенсивно развиваются.

В то же время появился новый класс сложных задач, к которому относятся расчёты крупных электрофизических установок, например, для физических экспериментов, готовящихся на Большом Адронном Коллайдере (LHC) в Европейском Центре Ядерных Исследований (CERN, г. Женева, Швейцария). Одной из таких установок является детектор ALICE [8]-[11]. Одновременно в проектах российских физических экспериментов [12, 13, 14] главным фактором становится минимизация затрат на создание установок, что стимулирует применение оригинальных конструктивных решений и, как следствие, интенсивное использование математического моделирования с развитием известных методов расчёта и созданием новых алгоритмов. При этом сложность моделирования электрофизических устройств часто связана с необходимостью учёта поведения векторных полей различной природы.

Можно выделить следующие принципиальные трудности, возникающие при решении подобных задач с помощью обычных подходов: поведение стационарных векторных полей описывается, как правило, системой дифференциальных уравнений с операторами векторного анализа, такими как дивергенция, ротор и градиент. Во многих случаях введение скалярного или векторного потенциалов приводит к сокращению числа уравнений системы. Однако при таком подходе необходимо следить за точностью последующего вычисления карты поля посредством дифференцирования потенциалов, поскольку операция численного дифференцирования, вообще говоря, является некорректной [15]. Ещё один недостаток такого подхода проявляется, когда поведение искомых векторных полей в разных областях задаётся разными уравнениями и на границах этих областей должны выполняться определенные условия для векторов поля. При использовании потенциалов о выполнении этих условий можно говорить в смысле сходимости приближённых решений к точному, то есть только в пределе; при решении задач относительно векторов поля, без введения потенциалов, обычно используются смешанные формулировки, после дискретизации которых возникают системы с так называемой седловой *«. точкой. Итерационные методы решения таких систем для сложных задач пока недостаточно эффективны; алгебраические системы со специфическими свойствами возникают при использовании некоторых алгоритмов для точного учёта убывании векторного поля на бесконечности и применении некоторых схем метода -Галёркина с разрывными базисными функциями. Для сложных задач такого типа также нет достаточно эффективных итерационных методов; очень важным для практических задач является вопрос о реальной точности вычисляемых полей при использовании ячеек конечных элементов с границами, которые задаются поверхностями второго порядка. Обычно точность в этом случае зависит от геометрических свойств ячейки, которыми определяются оценки для якобианов прямого и обратного преобразования в стандартный элемент. Поэтому использование базисных функций высокого порядка аппроксимации для потенциала в общем случае не гарантирует высокую точность аппроксимации векторного поля; и, наконец, довольно часто остаётся открытым вопрос о точности описания сложного электрофизического устройства уравнениями для проведения расчётов.

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

2. Общая теория проекционно-сеточных методов для эллиптических задач с монотонными операторами.

Разработка теории проекционно-сеточных методов восходит к работам

B.Ритца (1908г.), И.Г.Бубнова (1913) и Б.Г.Галёркина (1915г.), причём метод Бубнова-Галёркина (взвешенных невязок) обладает очень большой общностью. Отметим также книги по проекционно-сеточным методам

C.Г.Михлина, Х.Л.Смолицкого [16], М.А.Красносельского, Г.М.Вайникко и др. [17], Р.С1аг1е1 [18], Л.А.Оганесяна, Л.А.Руховца [19], Г.И.Марчука, В.И.Агошкова [20], В.В.Шайдурова [21], СЬ.Б^аЬ [22], O.C.Zienkiewicz, 11.Ь.Тау1ог [23], Р.Мопк [24], Р.ёоИп [25]. Общие подходы к решению нелинейных эллиптических задач были заложены Л.В.Канторовичем [26], М.М.Вайнбергом [27], Р.Е.В]^с1ег [28], [29], М.К.Гавуриным [30]. Отметим также работы Е.П.Жидкова, И.В.Пузынина [31], а также М.М.Карчевского, А.Д.Ляшко [32] и других представителей этих школ.

В настоящем разделе приводятся основные положения теории для нелинейных задач с монотонными операторами. Используемые далее обозначения соответствуют книгам [33, 34, 21], а также близки к принятым в книгах [27, 35, 36]. Следуя [33], введём следующие определения.

Определение 1. Оператор А £ (X —■> X*) называется: радиалъно непрерывным, если при любых фиксированных u,v £ X вещественная функция s —» -f si;),?;) непрерывна на [0,1]; деминепрерывным, если из ип и в X следует Аип —^ Аи в X*; липшиц-непрерывным, если существует такая постоянная М, что M\\u-v\\ для любых u,v £ X.

Определение 2. Пусть u,v - произвольные элементы из X. Оператор A £ (X —» X*) называется: монотонным, если (Au — Ау, и — v) > 0; строго монотонным, если — Av. и — v) > 0, и ф v\ сильно монотонным (с постоянной монотонности т), если

Au — Av, и — v) > т\\и — г»(|2, т > 0.

Определение 3. Оператор А £ (X —> X*) называется коэрцитивным, если существует определённая на [0, оо) вещественная функция 7 с lim 7(5) = +00, такая, что (Аи,и) > 7(|Ы|)|М|. s—>00

Определение 4. Говорят, что оператор А £ (X —X*) обладает (S) -свойством, если из ип и в X и (Аип — Аи, ип — и) —> 0 следует, что ип и в X.

Необходимо отметить, что для монотонных операторов А £ (X —> X*) понятия радиальной непрерывности и деминепрерывности совпадают [33].

Для монотонных операторов справедливы следующие теоремы существования.

Теорема 1. [33] Пусть А £ (X —■> X*) - радиально непрерывный монотонный коэрцитивный оператор. Тогда множество решений уравнения Аи = f при любом / £ X* непусто, слабо замкнуто и выпукло. Теорема 2. [33] Пусть оператор А £ (X —> X*) радиально непрерывен, строго монотонен и коэрцитивен. Тогда существует А~г £

X* —> X) и этот обратный оператор строго монотонен, ограничен и деминепрерывен. Если оператор А обладает, кроме того (S)— свойством, то А"1 непрерывен. Теорема 3. [34] Пусть

1) Н - гильбертово пространство со скалярным произведением (u,v),

2) Т - оператор из Н в Н,

3) существуют т > О, М > 0; М > т, такие, что

Ти — Tv\\h < М\\и — г>||я>

4) для всех u,v Е Н справедливо неравенство

Ти -Tv,u~v)> m\\u - v\\2H.

Тогда уравнение Ти = / имеет ровно одно решение для любого f Е Н и его можно получить методом последовательных приближений.

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

Лемма 1. [33] Пусть X - гильбертово пространство и оператор В G (X —> X) такой, что

Вх — By, х — у) > т\\х — у\\2, т> 0, \\Вх - Ву\\ < М\\х - у\\, для каждого х,у 6 X. Тогда оператор Ut Е (X —> X), определенный формулой

Utx = x- tBx, Ух Е X, является сжимающим при t Е (0,2т/М2) и

Utx - Uty\\ < Щ)\\х - у\\, где k{t) = (1 - 2mt + МЧ2)1'2 < 1.

Другие теоремы об итерационном решении нелинейных уравнений в гильбертовых пространствах, том числе для дифференцируемых операторов, можно найти в книге [37]. Следует также отметить работы [30] и [31], где в итерационную схему метода Ньютона и его модификаций вводится непрерывный параметр. Впоследствии они послужили основой для многочисленных публикаций с подобными итерационными схемами.

При решении практических задач важным понятием является потенциал оператора задачи, имеющий, как правило, определённый1 физический смысл. Приведём несколько теорем из [33, 34], которые использовались в диссертации для потенциальных операторов. Сначала напомним некоторые определения из [33].

Определение 5. Функционал F называется конечным, если F{x) € (—со, +оо) для каждого х е X.

Конечный функционал называется дифференцируемым по Гато, если существует такой оператор А 6 {X —> X*), что для любых х,у Е X lim ^{F{x + ty) - F{x)) = (Ах, у).

При этом А называют градиентом функционала F. Оператор А 6 {X —» X*) называется потенциальным, если существует такой конечный функционал F, что А является его градиентом. В этом случае F называют-потенциалом оператора А.

Определение 6. Функционал F называется выпуклым, если для любых х, у G X и t <Е [0,1] имеем tF(x) + (1 - t)F{y) > F(tx + (1 - t)y), и слабо полунепрерывным снизу, если F(x) < lim inf F(xk) для x G X и любой и—юоk>n последовательности {жп}^=1 в X, такой, что хп —х в X.

Справедлива следующая теорема о существовании минимума функционала. Теорема 4. [34] Пусть

1) X - рефлексивное, банаховое пространство;

2) F - функционал над X, такой, что lim F(u) = сю;

IMIx->оо

3) либо

I) F - непрерывный и выпуклый на X, либо

II) F - слабо полунепрерывный снизу на X. Тогда: a) inf F(u) > —оо; иеХ b) существует по меньшей мере один элемент щ £ X, такой, что

F(uq) = inf F(u)c) при выполнении условия пункта 2) и строгой выпуклости на X функционала F существует единственный элемент щ, обладающий свойством (Ь).

Приведём лемму о связи выпуклых функционалов и монотонных операторов.

Лемма 2. [33] Пусть конечный функционал F обладает градиентом А £ (X —> X*). Тогда следующие утверэюдения эквивалентны: a) функционал F выпукл; b) для любых i,j/6l функция s —> t/j(s) = F(x + sy) выпукла; c) оператор А монотонен; d) F(y) > F(x) + {Ах, у — х) для каждого х £ X и каждого у 6 X. Теорема 1 для потенциальных операторов формулируется в следующем виде.

Теорема 5. [33] Если А £ {X —> X*) - монотонный коэрцитивный потенциальный оператор, то уравнение Аи = / при любом f £ X* имеет решение.

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

Лемма 3. [33] Пусть А £ (X —> X*) - потенциальный оператор с потенциалом F. Для того чтобы элемент и £ X удовлетворял уравнению Аи = f, f £ X*, достаточно выполнения условия F{u) — и) = тт(.Р(г|) — (/, v)). Если оператор А вдобавок монотонен, то это vex условие также и необходимо.

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

Сначала определим понятие галёркинского решения, следуя [33]. Через {hn} С X обозначим полную систему линейно независимых элементов в X. п

Определение Т. Элемент ип G X вида ип = Y1 ^Г^г называется пг=i м решением Галёркина ( или галёркинским решением) уравнения Аи — / (относительно системы {hn}), если вектор {а",- --,0^} € Rn является решением, вообще говоря, нелинейной системы уравнений

Aun,hj) = (/, hj), j = l.,n.

Приведём основные теоремы о сходимости галёркинских решений. Теорема 6. [33] Пусть А б (X —> X*) - радиалъно непрерывный строго монотонный коэрцитивный оператор. Тогда при любом п существует точно одно решение Галёркина ип и ип —^ и в X. Это и является (единственным) решением уравнения Аи = /.

Теорема Т. [33] Пусть А € {X —> X*) - радиально непрерывный строго монотонный коэрцитивный оператор, обладающий (S) - свойством. Тогда при каждом п суш>ествует точно одно решение Галёркина ип и ип —* и в X. Это и является (единственным) решением уравнения Au = f.

Через Хп, следуя [33], обозначим линейное подпространство в X, натянутое на h\,., hn, с нормой, индуцированной из X. Для уравнений с сильно монотонными и липшиц-непрерывными операторами справедлива следующая теорема ( аналог леммы Cea для линейных эллиптических задач [18]).

Теорема 8. [33] Пусть A G (X —> X*) сильно монотонен и липшиц-непрерывен. Тогда, в обозначениях теоремы 7, имеет место оценка ип-и\\ <—d(Xn,u), d(Xn,u) = inf ||v —u||, m vexn где М - постоянная Липшица, а т - постоянная монотонности для А.

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

Отметим связь метода Галёркина и метода Рица.

Определение 8. Пусть А Е (X —> X*) - потенциальный оператор с потенциалом F. Элемент ип Е Хп называется п-м решением Рица для уравнения Аи = /, если F(un) - (f,un) = min^^) - (f,v)). veXn

Лемма 4. [33] Пусть A E (X X*) - монотонный оператор с потенциалом F. Элемент ип Е Хп точно тогда служит п-м приближённым решением Рица уравнения Аи = /, когда он является п-м галёркинским решением этого уравнения.

В пакетах и комплексах программ (например, [38, 39]), предназначенных для решения эллиптических задач, для дискретизации дифференциальных уравнений применяют в основном метод конечных элементов и реже -метод конечных разностей, а для дискретизации интегральных уравнений -метод коллокаций с кусочно-полиномиальной аппроксимацией неизвестной функции. Иногда используют сочетания конечных и неограниченных ("бесконечных") элементов, а также конечных элементов и граничных интегральных уравнений.

В диссертации рассматриваются только узловые скалярные и векторные конечные элементы. Информацию о рёберных и граневых (facet) конечных элементах можно найти, например, в [24, 25]. Следуя книгам [21, 18], приведём основные теоремы об интерполяции конечными элементами. Сначала определим необходимые понятия в соответствии с [21]. Определение 9. Конечный элемент в Rn есть тройка (w, Р, Е), где

1) w Е Rn - замкнутое подмножество с липшицевой границей и непустым множеством внутренних точек;

2) Р - N -мерное пространство функций, определенных на и>\

3) Е - набор линейно независимых линейных функционалов (рг '• Р —> г = 1, .ТУ.

При этом предполагается, что набор Е является Р - разрешимым в следующем смысле: для любого набора вещественных щ, г = существует единственная функция р Е Р, удовлетворяющая условиям у?г(р) = щ, Уг = 1, .ДГ.

Это равносильно существованию такого набора функций рг Е Р, г = 1,., ТУ, что

У 1)3 = 1, .А7".

С помощью этого базиса любая функция р £ Р представима в виде n

Р = ^2^г{р)Ргг=1

Линейные функционалы называются степенями свободы. Обычно </?г(р) - это значения функции р или её частной производной в некоторых точках а3 бш,]= 1, .ш, называемых узлами. Функции рг называются базисными функциями конечного элемента.

Часто конечным элементом называют само подмножество ги, поскольку с ним связаны как узлы для так и область определения

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

Определение 10. Интерполянтом функции V Е С3{уо) на конечном элементе (и>,Р, £) называют представление n

1и>У = У^<Рг(у)Рг, 2 = 1 где Рг - базисные функции в Р. Ввиду Р -разрешимости набора Е, это представление единственно.

Определение 11. Говорят, что для элемента (w, P. Е) выполнено условие полноты с целым к, если существует такое целое число к, что Р Э Р^(ги). то есть пространство Р допустимых функций содержит множество Pk{w) полных многочленов степени к, определенных на w.

Определение 12. Говорят, что для элемента (w, Р, S) выполнено условие регулярности, если он получен из элемента (и>,Р, Е) с помощью взаимно однозначного преобразования координат F : у х, х,у Е Rn, при котором:

1) w переходит в w = {F(y), у G w}]

2) Р переходит в Р = {р(х) — p(F~1(x)),p 6 Р};

3) Е переходит в Е = {<£г(р) = ^ — 1; •■■> " и существуют константы С\, .,С4, такие, что выполняются неравенства cihr\v\r>{S < \F~lv\uw < c2hr\\v\\r^,

О < С3 < sup J/ inf J < C4. w w

Здесь v € W^w), Л- - диаметр ги, a J - якобиан преобразования F. Теорема 9. [21] Пусть условие полноты для элемента (w, Р, Е) выполнено с целым к, таким, что вложено в Cs(w), а для элемента (uj,P, Е) выполнено условие регулярности для всех т < k + 1. Тогда для любого v W

Ik - IwV\\r,w < c(w, P)hk+1-r\\v\\k+1>{d.

При практических расчётах в качестве F обычно используют афинные и изопараметрические преобразования. Для афинных, более простых преобразований, оценки погрешности интерполяции значительно упрощаются [21, 18], и оценка теоремы 9, справедливая для норм, может выполняться и для полунорм тех же пространств. Для наиболее распространенных двумерных конечных элементов детальные оценки в случае афинных и изопараметрических преобразований получены в [20].

Оценка погрешности интерполяции для одного конечного элемента при определённых условиях переносится на всю расчётную область Г2. Для этого рассматривается совокупность конечных элементов (и>к, Рк, Ек), к = 1,.,К1г, покрывающих область О. В этом случае при естественных предположениях [21] об узлах общей сетки, набора линейно-независимых функционалов 3 = 1, и получающихся для каждого узла сетки базисных функций г — 1,., Л/71, справедлива следующая теорема. Теорема 10. [21] Пусть все конечные элементы (т^, Рк, Т,к), к = регулярны с показателем к + 1 > г с представителями некоторого множества элементов Фе = {(и»,Р, 2)}, все элементы которого удовлетворяют условию полноты с равномерно ограниченной нормой оператора 1Ш : И^^)- Тогда для любого V £ при условии вложения И7^4"1 в С5(0) справедлива оценка ; у-1Н^\\г,п < сНк+1~г\\у\\к+11П, где с - положительная константа.

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

3. Краткий обзор методов решения эллиптических задач с контролем точности

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

В первом подходе используются конечные элементы с убывающими значениями диаметров, но с фиксированной степенью многочленов, входящих в функции формы. Этот подход реализован в ряде программ, предназначенных для решения линейных и нелинейных эллиптических задач. К ним относятся, например, [40] (двумерный случай) и известное семейство программ ANSYS [41]. Содержащиеся в них адаптивные алгоритмы [42] позволяют пользователю в полуавтоматическом режиме получить приближённое решение с контролируемой точностью в большинстве практических задач.

Во втором подходе сгущение сетки производится за счёт наращивания степени многочленов функций форм. Здесь в качестве примеров можно привести программы, предназначенные для решения задач теории упругости: [43] (линейный двумерный случай), MSC/PROBE [44], PRO/MECHANICA [45] и некоторые другие. В них также используются адаптивные алгоритмы с полуавтоматическим режимом.

Третий подход объединяет два первых [46, 47]. Его основным достоинством считается экспоненциальная скорость сходимости приближённых решений к точному решению задачи (по крайней мере для некоторых модельных задач) [46, 48] в зависимости от числа степеней свободы. Примером программы, реализующей этот подход, служит двумерная программа HERMES [49], в которую включены адаптивные алгоритмы с иерархическими рёберными конечными элементами. Ещё одним примером является программа 3Dhp90 [50], предназначенная для решения трёхмерных задач. В ней используется специальная hp-адаптивная стратегия [51].

Все адаптивные алгоритмы основаны на генерации последовательности приближённых решений. Для получения каждого элемента данной последовательности естественно использовать самые быстрые методы решения систем линейных алгебраических уравнений. К ним можно отнести многосеточный (тиН^пс!) метод [52, 53, 54, 21], метод многоуровневых предобусловливателей [55], некоторые методы декомпозиции области [48]. При этом важным фактором, влияющим на выбор того или иного алгоритма, является возможность его эффективного распараллеливания и векторизации.

Для сравнения трех вышеперечисленных подходов приведём оценки скорости сходимости приближённых решений к точному решению задачи из обзора [46] для задачи Неймана в Л2:

Здесь П - многоугольник общего вида, / Е 1/2 (Г2). Для численного решения этой задачи используются различные семейства триангуляций Тд с некоторым минимальным углом во всех треугольниках, а также с фиксированными минимальными и максимальными диаметрами треугольников ктгП(Ть) и /гГггах(Т/г). Через 5(Г/1,р) обозначается множество функций из которые в каждом треугольнике семейства Т^ является многочленом степени р > 1. Далее для Ь-версии используется обозначение 5(Т/1,£>о), где ро - фиксированное число; для р-версии - 5'(Т/1о,р), где /¿о - фиксированный параметр. Предполагается, что имеет вершины Аг) г = 1,2,.,7п с внутренними углами а решение задачи может быть записано в виде У21А + 72 = /, жеП, ч ди/дп = 0, ж Е д£1.

777 г=1

Здесь щ Е И/22(^); для задания ф{ используется полярная система координат (г, в) с центром в каждой вершине А^ функция щ является урезанной частью бесконечно-дифференцируемой функции; оц = 7г/7^

Представим оценки скорости сходимости приближённых решений для двух множеств. Первое определяется следующим образом: т

32 = {ие \¥}(П) :и = и! + < Л XIс2 < В},

1 ¿=1 где А, В - некоторые константы. Второе имеет вид гк

1(3 = {и б #№), к > ||Ф0+|анДои|иа(П) ^ с ■ - О'.}. где (1 > 1,1 — 2, тп г=1 здесь п = (Ив^х.Аг), пространство Нкр1(0) состоит из элементов, ограниченных по следующей норме:

К|а|</с где, как обычно, = д^/дх^дх^2. Н = а\ + «2, а о;1,а2 - целые неотрицательные числа. Пусть при г = 1, 2

ЩЩТьр^Зг) = Ы вир Ы ||и - х\\]¥}(П)-Согласно [46], для Ь-версии справедливы следующие оценки: где 5 = тт{1,7г/7о}, 70 = тах(7г), а С^ро), Сг(ро) - положительные г константы, зависящие от ро. Если сетка достаточно хорошо раздроблена, то справедливы неравенства

Для р-версии имеем [46]

Я(А^(Т/го, р), £2) < С2 ) (ТУ , р)) - } 5 где постоянная (^(Т/^) зависит от диаметра элементов Соответственно для Ир-версии справедливы оценки [46]

ЩЩТьр))-1'2 < Я(^(ТЛ|р),Я2) < С2тп,р))-1'2,

ЩЩТърЪБЗ) < Сехр^-цУЩТьр)), где С - некоторая константа, постоянные С\,С2 зависят от Ат{Т}ир), а /1 > 0, но её величина точно не известна. В трёхмерном случае для рассмотреной задачи, когда 5з есть некоторое обобщение множества 63, для Ьр-версии имеет место следующая оценка сходимости приближённого решения [46]:

Д(ЛГ,5з) < Сехр(-/гУлГ).

Следует отметить, что оценки экспоненциального вида возможны только при выделении особенностей решения рассматриваемой задачи, если таковые имеются. В этом случае они справедливы и для р-версии метода конечных элементов. В [56, 57] были впервые исследованы решения двумерных эллиптических задач с угловыми особенностями. А экспоненциальные оценки, видимо впервые, получены в работах Е.А.Волкова [58, 59]. В качестве примера можно привести пакет программ ЭДИП [60], в котором предусмотрена возможность выделения угловых особенностей задачи при использовании метода конечных разностей.

В работе [47] дан типичный алгоритм адаптивного решения эллиптической задачи проекционно-сеточным методом. Он состоит из следующих шагов:

1) Ввести данные начальной сетки;

2) Решить дискретизованную задачу;

3) Оценить ошибку приближённого решения;

4) Если ошибка достаточно мала, то перейти к обработке полученного приближённого решения и закончить процесс решения. Если же ошибка недостаточно мала, то изменить сетку с целью уменьшения ошибки и перейти к выполнению шага 2).

Укажем некоторые особенности автоматического адаптивного hp-алгоритма для решения двумерных задач теории упругости, который описан в работе [43]. Этот более сложный алгоритм содержит такие шаги, как определение способа улучшения конечно-элементной сетки (имеется ввиду h- или р-улучшение), необходимой степени р для функций формы с целью достижения требуемой точности, а также необходимого дробления элемента вне окрестности точки с особенностью и в её окрестности. В этой же работе приведены формулы для оценок ошибки приближённого решения на конечном элементе, в том числе в окрестности точки с особенностью, а также даются необходимые ссылки на публикации, в которых эти формулы получены. Более современные алгоритмы рассматриваемого типа описаны в публикациях [51, 61, 62].

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

Рассмотрим задачу нахождения во всём пространстве распределения магнитного поля, созданного стационарными токами в проводниках и намагниченностью находящихся в пространстве изотропных ферромагнетиков. В случае отсутствия поверхностных токов и токов, протекающих по ферромагнетику, эта задача сводится к нахождению вектор-функций В, Н из системы уравнений Максвелла для стационарного магнитного поля (например, [63]):

V • Б = О, УхЯ-J, ё = (20 цН, (1) lim Н(х) = 0. (2) х—>оо

Здесь используются следующие обозначения: В, Н - вектора индукции и напряжённости магнитного поля; 3 - известный вектор плотности тока, отличный от нуля в ограниченной области Г^; ц = /л(\Н\) или ¡1 = ¡¿(\В\) - функция магнитной проницаемости ферромагнетика, задаваемая в ограниченной области Цр, вне этой области ¡л = 1; //о - магнитная проницаемость вакуума; х- точка пространства. На границе Г раздела сред с различными магнитными характеристиками выполняются условия: п-В] = 0, [п х Н] = 0, (3) где п - вектор нормали, а выражения [п • В] и [п х Н] означают разность пределов стоящих в скобках функций при стремлении точки наблюдения на границу Г изнутри и извне области Пр.

Далее для простоты изложения будем рассматривать случай двух сред - магнитной и немагнитной. Относительно свойств функции ¡1 обычно предполагают [64] выполнение неравенств

О < м* < М < АЛ 0<^< (1\В\/(1\Н\ < /Л (4) где - положительные константы.

Для численного решения задачу (1)-(3) обычно формулируют [38, 39] либо относительно векторного потенциала Л, который вводится по формуле В = V х А при условии V • А = 0, либо относительно скалярных —* потенциалов (р) вводимых по формулам Н = — \7у> в ^ и Н = Чч> + я5 вне области Пр. Через Н5 обозначают решение задачи (1)-(3) при отсутствии магнитных материалов (закон Био-Савара) где \х — у\ - расстояние между точками хну.

Задача (1)-(3) может быть также сформулирована относительно модифицированного или обобщённого скалярного потенциала [65].

Для векторного потенциала получают следующую задачу [66, 67]:

V х х I) - • А) = I,

- с условиями на границе раздела сред п х х А] = 0, [п • А] = О, п х А] ~ 0, [1/(^0^) V ■ А] =0. В случае двух скалярных потенциалов имеем задачу [68]

V • /¿(|У^|)У^> = 0, х € Пр; V • У(/? = 0, же113\П^ с условиями на границе раздела сред п • (м(¡V'0|)V-0 - + Я5) = 0, а также рх ф = + (5) у

Отметим, что при выполнии хотя бы одного из условий <фбГ = 0, или / ^сЮ = 0,

Уг где Г, Ор - соответственно части границ Г и уравнения (5) заменяются на условие п х {у-ф - Ус^ + Я5) = 0.

В работе [69] показано, что последнее условие предпочтительнее уравнений (5).

Задачу (1)-(3) можно сформулировать в интегральных уравнениях относительно векторов поля

Я - -^У / М ■ у2—¿П = Я50), 4тг \х-у | где М = Я//7,0 - Я.

Отметим, что условие (2) в дифференциальных формулировках обычно учитывается приближённо, а в интегральных уравнениях - точно. Поэтому для точного учёта условия (2) в дифференциальных формулировках к ним добавляется интегральное условие [70] на границе Гх, содержащей область

1 (* ^В п п(В ■ п) + п х (Я х п) + —V / 1—

2-л \х - у| ^-с1Ту = 2Н5{х). 2тг уГ1 \х-у\

С этой же целью в дифференциальную формулировку задачи (1)-(3) относительно потенциалов вводят граничное интегральное уравнение, например, для гармонической функции и и гладкой границы Г\ в виде

2тги(х) — [ -—-—гп • \7udSy — [ ип • V-.—-—^Бу,

7 \х-у\ ] \х-у\

Г1 г:

X 6 Гь

Это уравнение можно рассматривать как обобщение третьего краевого условия для эллиптических задач, которое можно включать в общую формулировку задачи и как граничное условие Дирихле, и как граничное условие Неймана [71].

Каждая из приведённых формулировок задачи (1)-(3) имеет свои достоинства и недостатки, и этим определяется её область применимости.

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

Как уже отмечалось, для решения задачи (1)-(3) необязательно находить потенциалы. Так, при решении интегральных уравнений относительно вектора В [72, 73, 74] методом коллокаций с кусочно-полиномиальной аппроксимацией неизвестных этот вектор находится в области Qp, а затем с помощью пересчёта он определяется в области Я3 \ Q.F- В работах [75, 76] рассматривается возможность решения задачи относительно векторов Н, В с помощью смешанного метода [21] с использованием "рёберных элементов "(edge-elements), обобщением которых являются элементы Неделека (Nedelec's elements) [77]. В работах [78, 79, 80] также рассматриваются формулировки задачи относительно векторов индукции или напряжённости магнитного поля. Причём, несмотря на то, что в подходе, описанном в [79], задача решается относительно вектор-функции, конечно-элементная схема сохраняет хорошие свойства соответствующих схем для потенциалов.

Доказательству разрешимости вышеописанных формулировок, которое обычно проводится с использованием неравенств (4) и теории монотонных операторов-, посвящены работы [81, 38], [82, 83, 84], [70, 85, 86]. В [38] приводится обширная библиография до 1990 года. Результаты по общей теории смешанного метода содержатся в [21], а также, например, в [87, 88].

Поскольку задачу (1)-(3) приходится решать при проектировании практически любого ускорителя, а также при разработке спектрометрических электромагнитных систем для физических экспериментов, при расчётах электромоторов, трансформаторов, громкоговорителей и другой бытовой техники, то естественным является появление различных пакетов и комплексов прикладных программ, предназначенных для её решения. Например, в обзоре [38] и в книге [39] представлен широкий перечень программ, которые используются в физических институтах России и зарубежом. В работе [89] дано описание комплекса программ с учётом применения векторных и параллельно-векторных компьютеров. За прошедшие 10-15 лет появились новые программы (например, [90. 91]), однако, как отмечено в [94], существующий уровень стандартов информационной поддержки наукоёмких изделий в промышленном производстве развитых стран мира пока остаётся таким же, как и в последнем десятилетии прошлого века.

5. Краткое содержание диссертации

Первая глава посвящена проекционно-сеточным методам решения краевых задач для векторов с интегральными граничными условиями, точно учитывающими убывание решения на бесконечности. Как известно из работ Р.П.Федорепко [52] и Н.С.Бахвалова [53], особенностью численного решения линейных эллиптических задач в дифференциальном подходе является быстрая сходимость метода простой итерации для быстро осциллирующих компонент решения задачи. В случае же интегрального подхода, как показано, например, в работе Е.П.Жидкова, Б.Н.Хоромского, О.И.Юлдашева [95], в методе простой итерации быстрее сходятся гладкие компоненты решения задачи. В связи с этим комбинированный подход, который объединяет дифференциальный и интегральный подходы и позволяет точно учитывать условие убывания решения на бесконечности, продолжает вызывать интерес. Среди первых публикаций, посвящённых разработке комбинированных методов, можно отметить работы [96, 97, 98, 99, 100, 101, 70, 102]. По-прежнему разрабатываются прикладные пакеты и программы на эту тему [103, 104], что свидетельствует об актуальности излагаемых здесь результатов.

В §1.1 доказывается теорема 1.1.1 о компонентах непрерывно-дифференцируемого вектора. Как известно, каждый вектор можно разложить по направлению единичного вектора и ему перпендикулярному. В силу теоремы Гельмгольца [105] единичный вектор может быть градиентом, ротором или их суммой. Поэтому при описании компонент вектора в новом разложении получаем различные системы уравнений для его определения.

Во-первых, из теоремы 1.1.1 вытекает абстрактная формулировка внутренней краевой задачи

V ■ (<31 (х, й)и) = /(х, и), х Е

V х (а2(х,и)и) = Р(х,и), х Е Г2; (6) а{х)п(а\й • п) + (3(х)п х [а2й х п) = д\, х Е дП, (7) где все функциональные зависимости предполагаются заданными, а во-вторых, задача вида

V • (А\(х,и) х и) = ж е

У(А2(х,и) ■ и) + V х (Л3(ж,й) х и) = Р{х,й), ж 6 П; (8) а(х)п{(А1 хи)-п) + [А2 ■ и)) + (3{х)п х (А3 х и) х п) = д2, хЕ <90. (9) Также из теоремы 1.1.1 следует формулировка задачи

V • (ах(ж, й)и) + V ■ (А\(х, и) х и) = /(ж, и), х Е О;

Vx(a2(x}й)й)+V(A2(x,й)^u)+Vx(A3(x,u)xй) = Р(х)и), х Е П; (10) а(а:)п((ахи • п) + (А\ х и) ■ п) + (А2 ■ и)) + (3{х)п х ((а2й + х г?) х п) = (11)

Все три абстрактные задачи (6),(7), (8),(9) и (10),(11) при определённых коэффициентах и правых частях имеют смысл, поскольку уравнения каждой задачи описывают поведение компонент непрерывно-дифференцируемого вектора.

Сформулируем абстрактные задачи для описания стационарных векторных полей во всем пространстве. Если дифференциальные уравнения (6) дополнить условиями

ИтЬЗ = 0, г = 1,2, (12) ж|-»оо |ж| то имеет смысл задача (6),(12). А при условии на бесконечности

Hm üiiiil = о, « = 1,3, = 0 (13) х\—>оо кГ |х|—>оо имеет смысл задача (8),(13). Аналогично, имеет смысл задача во всем пространстве (10),(12),(13). Для краткости дальнейшего изложения рассматривается частный случай задачи (6),(12) в следующем виде:

V • (а\(х,и)й) — f(x), V х (а2{х,й)и) = F(x), х Е П;

V- (oi(a;)it) = Л(ж), V х (äiWu) = Fi(x), x6fii = (R3\^); (14) n • (aiü — ä\u) =0, n x (CL2Ü — ä\u) = 0, x E Г = dQ; lim а\й = 0. x|—»oo

Здесь граница Г, разделяющая области Q и f2i, является липшицевой, / Е L2(ii), fi Е L2(i2i), F Е (Ь2(П))3, Fi E i))3, все коэффициенты уравнений предполагаются непрерывно-дифференцируемыми в Q. U

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

Частным случаем задачи (14) являются задачи магнитостатики. Здесь следует отметить, что в автореферате кандидатской диссертации автора [71] для двумерных задач магнитостатики было предложено рассматривать граничные интегральные уравнения как специальные граничные условия (лемма 3), и на основании лемм 4,5 были обоснованы два метода включения граничных интегральных уравнений для векторного потенциала в конечно-разностную сеточную схему решения дифференциальной нелинейной задачи магнитостатики относительно векторного потенциала. По всей видимости, независимо от этих исследований, один из упомянутых методов включения, только в конечно-элементную схему, был реализован в виде специальной программы для решения подобного класса задач [103]. Впоследствии она использовалась для расчётов магнитного поля главного дипольного магнита Большого Адронного Коллайдера [104].

Следует отметить, что излагаемый в этой главе подход (§ 1.4) является более общим, чем методы, изложенные в работах [98, 99, 102], а также более экономичным по затратам памяти ЭВМ, чем подходы из [100, 101, 103].

Для исследования системы (14) удобно ввести обозначения у= . —' и = у-г. х Е fii, а\й, х Е Oi

Для х Е Q будем использовать следующие зависимости: Z = v(x,V)V и V — к(х, Z)Z.

Для решения системы (14) при / = 0 рассмотрим интегральное уравнение относительно векторного потенциала А, который вводится по —* —* формуле V = V х А.

Векторный потенциал удовлетворяет уравнению

RA = A — J- J U (А) х V^dy = Af{x), х Е П. (15) п

Здесь U(А) = (1 — и(х, V х А))\7 х А, а г = \х — у| - расстояние между точками х и у, Af(x) - заданный вектор, Af Е (Lг(П)) . Предполагается, что 0 < v* < v < 1, — const.

В § 1.2 исследуется интегральное уравнение (15) для векторного потенциала в гильбертовом пространстве

H(rot, tt) = {u:uE (Ь2(П))\ V хиЕ {Ь2{П))3} со скалярным произведением и, v)a = at(u, v) + (V х й, Vxi), где а — const, а > 0, и, v 6 H(rot,fl), (•, •) - скалярное произведение в

Через || • || и || • || а обозначим соответствующие нормы в пространствах (L2(^))3 и H(rot,Q).

Основным условием, которому должны удовлетворять коэффициенты системы (14), является следующее неравенство:

Ui - Ü2\ < (1-^)|У1 -V2\, хеп. (16)

Свойства оператора R : H(rot, О) —H{rot, Г2) устанавливаются в следующей теореме. Теорема 1.2.1

Если выполняется условие (16), то:

1) R - сильно монотонный оператор, то есть для Äi,Ä2 £ H(rot,Q)

Л(^) - R(Ä2),Ä1 - Ä2)a > XI IIA - Ä2\\l, где а, xi ~ положительные константы;

2) R удовлетворяет условию Липшица, то есть для А\,А2 £ H(rot, П)

ЦЩЛг) - R(Ä2)\\a < kWh - Ä2\\a, где а, 1\ - положительные константы.

В теореме 1.2.2 устанавливаются свойства производной Фреше оператора R:

R'(Ä)Ä+ = Ä+- QU', где Ü' = (Е — Z'(V))V+, Z'(V) = {dZt/dV3}lJS=1,2fii V+ = V x Ä+. Формулируются условия при которых решение уравнения (15) можно получить с помощью итерационного процесса Ньютона.

Для дискретизации интегрального уравнения используется метод Галёркина. Введём пространство

Jl(ü) — {и ■. й £ (И^П))3, V • и = 0, й • п\г = 0}, где и • п|г - след на границе Г, a W^iP) - пространство Соболева со скалярным произведением й, v)ji = (V х и, V х v) и нормой Undji = ||V х гГ||, где и, v Е 1/1(Г2).

В леммах 1.2.1 и 1.2.2 доказываются вспомогательные результаты об эквивалентности норм и свойств оператора V х для элементов пространства Jx(fi). Далее рассматривается уравнение

PAR(A) - Af) = О, и доказывается теорема 1.2.3. Теорема 1.2.3

Если выполняется неравенство (16), то:

1) PjiR : —> «71 (Г2) - сильно монотонный и липшиц-непрерывный —* 1 оператор, то есть для А\, А2 Е J (О) выполняется неравенство (Pji(R(Ai) - R(A2)), Аг - A2)ji > - 12||2j1; а также

Рл(Я(Л) - < (2 - v*)\\A\ - l2||Ji;

2) в предположениях теоремы 1.2.2 оператор PjiR'(A) : J1 (17) —► Jx(f2) непрерывно обратим и удовлетворяет условию Липшица по А. При этом для А, А+ £ Jx(f2) выполняются следующие неравенства:

Pj,R!(A)A+\\j,>^\\A+\\j4

РАР!{Аг) - R'(A2))\\ji < 7|Hi - A2\\ji, где 7 = const.

Согласно методу Галёркина, уравнению (15) соответствует система

PJL(R(Am) - Af) = О, где J^ - m-мерное линейное подпространство в J1(i2) с нормой, индуцированной из t71(Q), базисом щ, .,nm, и

Лд тп агиг. г=1

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

Из определения нормы J1 следует, что уравнение (15) эквивалентно удобной для численной реализации комбинированной m-мерной системе

J uV х Am • V х utdx + J n • V x иг J ~—(Ujm • n)dSydSx — fl г г J(n x щ) ■ V x AfdS, i = 1, 2, .ш; г —t —* здесь Ujm = (1 - Pjmv)V x Am, а проекция PjmvV x Am находится из специального уравнения. Jm - конечномерное подпространство в пространстве J = {и : и Е V • и = 0}.

В § 1.3 исследуется интегральное уравнение для скалярного потенциала. —*

В случае, когда F = 0, в задаче (14) можно ввести скалярный потенциал if по формуле Z — V</9. Он удовлетворяет уравнению

T(ip) = ф) + j-J U((/?) . Vl-dy = iff(x), хеп. (17) п

Здесь U(f) = (k(x,Vip) — l)V</>, a iff - задаршая функция, iff б ^(П). Предполагается, что 1 < с к j /с * = const. Пусть C?i(n) - гильбертово пространство: Gi(fl) = {и : и = V^,^ £ W^Q)}, через Aq : G\(Q,) —» Gi(fl) обозначим интегральный оператор

А0и)(х) = V J u-V^dy. п

Известно, что оператор Ао самосопряжён и положительно определён, кроме того,

Ai < Що|| < Л2 < 1, Лк = const, > 0, к = 1, 2.

Учитывая эти свойства оператора Ао, в Gi(Q) удобно ввести скалярное произведение п,г7)* = (A^u, г>). Для соответствующей нормы имеем

HI < ми < ^И

В пространстве W21(H) ввёдем скалярное произведение с/?, VO^ =/%, + (V<A W)*, где (3 — const, (5 > 0, с/?, ф G И^П). Соответствующая норма || ■ ||| = (•, •)/?. Далее предположим, что выполняются неравенства

U1-U2W < (к* - l)\\Zi - Z2\\, (18)

Ui - U2, Z! - Z2) >0, (19) где Ui, Z{ & (Z/2(i2))3, ¿ = 1,2.

Справедлива следующая теорема. Теорема 1.3.1 Если справедливы неравенства (18) и (19), то:

1) Т является сильно монотонным оператором, то есть для любых

Т(ср 1) - T(<p2),ip 1 - 4*2)0 > хЛч>\ - Ч>2\\% где Р, Х2 - положительные константы;

2) Т является липшиц-непрерывным оператором, то есть для любых

T(ipi) — Т((/?2)||/з < h\\fi ~ V2II/3, где ¡3,12 - положительные константы.

Для решения уравнения (17) методом Ньютона определяется производная Фреше

Г(<р)<р+ = + <у\г) -1 п где = ^'(2") = В теореме 1.3.2 формулируются условия, при которых решение уравнения (17) можно получить с помощью итерационного процесса Ньютона.

Для дискретизации уравнения (17) по методу Галёркина рассматривается гильбертово пространство Я"1 (£7), состоящее из элементов пространства И^П), для которых выполняется хотя бы одно из условий

J ф(1Б = 0, или J ф<11 = О, Го ¿о где Го - часть границы Г, а - незамкнутая кривая в области П, со скалярным произведением

Пусть Нт - т-мерное линейное подпространство Н1 с базисом ф\,., фт и нормой из Н1. Тогда для аппроксимации Галёркина т г=1 получим уравнение

Рн 1 - <?/) = 0, (20) т где через Р^ обозначен оператор проектирования Н1 на

-"т

Уравнение (20) эквивалентно ш-мерной объединённой системе

J кЧфт-Уф1(1х + J пхЧфг J -^(и{с1)т><п)(13уС13х = п г г J'фin■'V(pfdS, г = 1,.,т, г где ЩС1)т = {Р(С1)тк - 1 а проекция Р(С1)т/сУ(/?гп находится из специального уравнения.

Учитывая свойства операторов проектирования, доказываются теоремы 1.3.3 и 1.3.4.

Теорема 1.3.3

При выполнении условий теоремы 1.3.1 операторы Р^Т : Н1 —> Н1 и Р^ 1Т : Н^ —»• являются сильно монотонными и липшиц-непрерывными. Теорема 1.3.4

При выполнении условий теоремы 1.3.2 производные Фреше Р^\Т' : Н1 Н1 и Ртгг Т' : Н1 —> Н^ являются положительно определёнными п т и липшиц-непрерывными.

Из теорем 1.3.3 и 1.3.4 следует не только существование и сходимость приближённых решений к точному, но также сходимость метода последовательных приближений и метода Ньютона.

В § 1.4 рассматривается класс комбинированных систем уравнений для нелинейной задачи (14), полученный с помощью введения вспомогательного вектора. Для простоты изложения предполагается, что граница Г является более гладкой и в каждой точке имеет непрерывный вектор нормали. —>

Для векторов V и Z имеем

V х ^ = у-У = /, У = хеп, (21) где функция к задана либо как к = к(х,У), либо как к = к(: —*

Используя формулу Грина в области для векторов V, Z, получаем д(х) = -V | + V х у ^ЛЗу + V, хе ПЬ (22) г г где п - вектор внешней к Î7 нормали. Здесь = v х /¿^ - v / t/y

Qi Oi

Учитывая скачки интегралов в (22) при ж-+Ги условия связи нормальной —* компоненты вектора V и тангенциальной компоненты вектора Z при переходе через границу Г, приходим к граничному соотношению

V ■ п fZ х п 1 ^ d^+Vx / dSy+-{n{V-n)+nx{Zxn)): х е Г. г г

23)

Здесь g1(x)=g(x)-±(nAf + AFxn) + V J ^-dSy + V х J П Х dSy, г г а Д/= /-/i, ДР = F — Fi.

Таким образом, в общем виде задачу (14) можно сформулировать как задачу с уравнениями (21) и граничным условием (23).

С целью получения в (23) только одного интеграла по Г вводится вспомогательный вектор-функция Р, удовлетворяющий в Г2 следующим уравнениям:

V х Р = О, V-P = 0, хеП. (24)

Тогда из формулы Грина имеем а(х)Р(х) = V х J ^^-dSy + V J ^^dSy, (25) г г где си(х) = 47Г, если х Е Г2, и а(х) = 2тт, если ж G Г. С учетом (18) и (20), задавая для вектора Р условия на Г в виде Р хп = Z хп или Р-п = V - п, получим следующие краевые условия к уравнениям (21) и (24):

Pxn = Zxn, (Р + V ^ ' ПdS) хп = д1(х)х n, (26) г Г а также г

Р • п — V • п, (Р + V х

О?-!!^) • п = дг(х) • п, (29) г точно учитывающие поведение решения на бесконечности. Каждая пара краевых условий (26)-(29) содержит только один граничный интегральный оператор и с уравнениями (21), (24) соответствует четырём комбинированным формулировкам.

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

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

Как известно, общие подходы к решению нелинейных эллиптических задач были заложены Л.В.Канторовичем, М.М.Вайнбергом, Р.Е.Вго"ш1ег, О.Л.М1г^у, М.К.Гавуриным и другими. Следует также отметить работы Е.П.Жидкова, И.В.Пузынина, а также М.М.Карчевского, А.Д.Ляшко и других представителей этих школ.

Можно выделить два направления в разработке и применении проекционно-сеточных методов с конечными элементами, которые особенно интенсивно развиваются в последнее время. К первому направлению относятся к—, р— и Кр— схемы (р > 2) метода Галёркина с разрывными базисными функциями [4, 106], а ко второму - h—,p— и hp— версии схем с векторными конечными элементами ( Whitny's, Nedelec's elements).

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

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

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

Z{Q) = {йе (Ж^П))3 : V • и = 0, V х и = 0}; и (О) = {йе (L2{Q))3 : V • и = 0, V х и £ (L2(ft))3}; = {гГ е (L2(f2))3 : V • и б L2(fi), V х и = 0}.

Функции из пространства Z(Q) могут быть выражены в виде градиентов гармонических функций. Поэтому в п. 2.1.1 приводятся алгоритмы генерации гармонических конечно-элементных базисных функций.

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

Н^г^М, (30) г—1 где иг, г = 1, .,т, - значения приближённого решения в узлах {жг}г=1,.,т на границе ди>, а А^(ж), г,= 1, .,т, - базисные функции.

Для определения гармонической базисной функции ^(х), то есть такой, что V • У7Уг- = 0, согласно первому алгоритму, используется метод коллокации. Она ищется в виде Щ(х) = а) /зСЖ^О) гДе " 1 неизвестные коэффициенты, д(]) -индексная функция, ^ = 1 , .,га, а в качестве ¡д(з) выбираются функции

Рп+1,к+1(Х) = ¿<пк(г/го)П С0Б(ку)Рк{с08б) ,

Чп+1,к+\{х) = (1пк(г/г0)път{ку)Рк(созв) ' (31) при х = (г, <9, <р), ¿пк — (2п + 1)(п — к)\/(п + &)!, которые вычисляются по реккурентным формулам и входят в общее представление решения задачи Дирихле для оператора Лапласа внутри сферы радиуса го- Неизвестные коэффициенты находятся в результате решения системы

ТП

1д(ЛЫ) = 6(Х^Х1)) I = 1, .,771. 1

Индексная функция д(э) подбирается так, чтобы система была разрешима. Точность аппроксимаций с помощью Щ(х) определяется степенями гармонических многочленов, которые входят в функции з = 1, .,т, и для которых аппроксимационная формула (30) точна.

По второму алгоритму Иг(х) = 1з(х)> ГДе fз выбираются

3=1 последовательно из функций (31), а неизвестные коэффициенты находятся в результате решения системы где Ьг - лагранжевая граничная базисная функция. Точность приближений функциями ЛЦя) зависит от точности аппроксимаций функциями Ьг, г = 1

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

В п. 2.1.3 представлены два алгоритма построения конечно-элементных векторных узловых базисных функций для аппроксимаций высокого порядка из пространств 2(0,), 17 (П). В обоих случаях используется тот факт, что как и сами гармонические функции, их градиенты обладают хорошими аппроксимационными свойствами.

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

3 т 3 т ад = = ЕЕ^Й^)' (32) к—1 з=1 к= 1 з=1 где ик,з - значения приближённого решения в узлах {х3}3=1}„ш}ТП на границе дш.

Согласно первому алгоритму, базисные функции ищем в виде а^'1^где а^ - неизвестные коэффициенты, з=1 гармонические функции из (31), д- индексная функция. Неизвестные коэффициенты находятся в результате решения системы при Х[ Е дш: Зга а г 3 ш д г Зт д л к^к1, кфк 2, к1 к2, 1 < к,к!,к2 <3, 1 = 1,2,.,т. (33)

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

Во втором алгоритме базисные функции ищются в виде ]¥^к\х) =

Зт . .

X] Щ уг где по правилу построения среднеквадратичного з=1 приближения коэффициенты Ь^'^ находятся из системы

3 7П р р г\ п

• V/í,(^)dS - / Ьг-^-аБ, 1 = 1,2, .,3т. (34)

О1 «/ */ ^

1- 5а; дш

Здесь г/^') = ^ + 1,^ = 1,2,., Зт, если = 1. Точность аппроксимаций в этом случае зависит от точности приближений лагранжевыми граничными базисными функциями Ь^ (1 < ъ < т).

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

Далее приводятся алгоритмы построения конечно-элементного базиса из гильбертовых пространств С/ (Г2) и У (Г2). Сначала рассматривается случай пространства У(Г2). Пусть приближённое решение и ищется в виде (32), где и^ - значения приближённого решения в узлах внутри и на границе ал

Согласно первому алгоритму, базисные функции представим в виде а^'^УНд(э){х), где а^ - неизвестные коэффициенты, д{э) з=1

- индексная функция, э = 1, .,т, - функции из множества

1; XI, х2,ж3; Х\Х2] хххъ\х2х3; хгх2х3-, х\, х\, х\,.}. (35)

Неизвестные коэффициенты находятся в результате решения системы, аналогичной (33). Во втором алгоритме базисные функции представляются в виде \¥^к\х) = ^ где по правилу построения 1 к I) среднеквадратичного приближения коэффициенты Щ ' находятся из системы, аналогичной (34). В обоих алгоритмах точность аппроксимаций базисными функциями зависит от степени многочленов (35), градиенты которых участвуют в определении коэффициентов а^'1^ и Ь^к,г\

Базисные функции из пространства II (£1) формируются из базисных функций Щк) 6 У(П), г = 1,., т, к = 1,2,3.

В п. 2.1.4 дано более подробное описание конечных элементов с гармоническими функциями формы из п. 2.1.1. В таблице 1 приведены такие характеристики конечных элементов, как вид ячейки ы, число узлов т и наибольший порядок с/ гармонических многочленов, которые точно приближаются соответствующими базисными функциями. Для стандартного симплекса 5 = {х1 > 0, г = 1,2,3,21 + 22 + 23 < 1} и пятигранной призмы Т = {21,22 > 0,21 + 22 < 1; —1 < 23 < 1} со стандартным расположением узлов на границах ячеек оба алгоритма дают одинаковые результаты. Однако, как видно из таблицы 2, в случае стандартного куба [—1,1]3 первый алгоритм демонстрирует более высокую точность.

Таблица 1 Таблица 2

5 Т т 10 20 34 52 18 38

1 2 3 4 5 2 4 алгоритм 1 алгоритм 2 т 26 56 98 26 56 98

3 5 7 2 3 4

Приведём пример сравнения сходимости приближённых решений в трёхмерной задаче интерполяции магнитного поля, когда в качестве конечных элементов используются последовательности обычных лагранжевых (параллелепипеды с 27, 64, 125 узлами) и гармонических (параллелепипеды с 26,56,98 узлами) элементов. На рис. 1 представлено поведение максимума модуля относительной погрешности с^1) в зависимости от общего числа узлов N, используемого при интерполяции на адаптивной сетке.

Как видно из рис.1, лучшая точность « 0.1 • Ю-4 достигается на 98 узловых гармонических элементах при общем числе точек интерполяции 16358. Это на 37 % точек меньше, чем при интерполяции 125 узловыми лагранжевыми элементами с меньшей точностью.

Рис. 1. Сходимость интерпо-7 ляций для параллелепипедов слева - направо, сверху - вниз): о - 27, 64, 125 узловые лагранжевые;

• - 26, 56, 98 узловые гармонические (первый метод).

12 О 5000 10000 15000 20000 2ЭООО N

В § 2.2 представлены обобщённые формулировки для численного решения краевых задач с системой из двух уравнений первого' порядка в дивергентной и вихревой формах относительно вектор-функций. Рассматриваются два частных случая задачи (6),(7), которые являются наиболее распространёнными на практике:

V ■ (аггТ) = /, = х е О; гГ = 0, хедП.

V ■ ш = /, Ух (а2й}) = Р, х е П; ^ ги = 0, х € дП.

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

Для построения обобщённой формулировки задачи (36) используется вариационный метод. Пусть а\ = 0,1(2:, у) - заданная функция, / 6 Ьг(О), юд(^)

Р Е (Ь2{^1)3 и V • Р = 0. Решение задачи ищется в пространстве И^о^)3: {и е (Ь2{П))3 : ди/дхк Е (Ь2{П))3: к = 1,2,3; й\да = 0}.

Как известно, такое пространство позволяет использовать две эквивалентные нормы: з

11^11^,0(0)3 = ^2{дй/дхк,дй/дхк)(Ь2тз-к=1

И1и?10(П)» = • + Х Й11(£Я(П))3

Очевидно, что решение задачи (36) обеспечивает минимум функционала хИ = (1/2) I[(V ■ (Й10) - Л2 + (V х у-Р)2}с1П. п

Предположим, что для коэффициента уравнения выполняются следующие условия:

С1|| V • < ■ (Й1гО|к(П) < ■ йЦ^п), (38)

V • (Й1(ж, ь^щ) - V • {ах(ж, 2)|и2(п) < с3\\щ - гУжу^з, (39) где сг > 0, сг = сопяЬ. г = 1, 2. 3, а V, щ, и2 £ Справедлива теорема. Теорема 2.2.1

Пусть = Ы)2 + (1 - í Е (0,1) и щ,у2 Е И^о^)3- Если для коэффициента ах выполняются неравенства (38),(39) и справедливо представление аг(х, -и+)у+ — ¿ах(ж, У2)У2 + (1 - £)ах(х, то при / Е 1/2 (П), Р £ (¿2(0))3 функционал определённый на И^Д^)3, является непрерывным, коэрцитивным и строго выпуклым. При выполнении (38), а также неравенства и (ах, у, и)г?)||^2(П) < с^И^уэдз, где u'(a\,v,u) = {д(аг(x,v)vl)/duj}itj=ii2,3) v,u G И^ДО)3, c4 = const, доказывается лемма 2.2.1 о том, что градиент функционала Fi, определяемый равенством n задаёт отображение из W^O)3 в (W^O)3)*.

При выполнении условий теоремы 2.2.1 и леммы 2.2.1 из известных теорем о связи метода Рица и метода Галёркина следует, что обобщённая формулировка для задачи (36) в виде уравнения

А^й) = о, Wg и/2уп)3, имеет единственное решение v G Ж^О)3- Кроме того, при каждом п это уравнение для Vun G (VF21io)n(^)3 имеет единственное галёркинское приближение йп G (И/211о)п(0)3.

Если в задаче (36) F — 0, то для её эффективного решения в полученной обобщённой формулировке рекомендуется использовать базис из пространства У (О), удовлетворяющий краевым условиям.

Далее вариационный метод используется для построения обобщённой формулировки задачи (37). Предполагается, что а2 = a2(x,w) - заданная функция. / Е 1/2(0), F G (L2(0))3 и V • F = 0. Решение задачи ищется в пространстве ^^(О)3. Функционал для задачи (37) имеет вид

F2(v) = (1/2) J[(V ■ (w) - /)2 + (V х (a2w) - F)2}dCt. n

Пусть для коэффициента уравнения выполняются следующие условия:

Csifv X й?||(£2(п))3 < ||V X (52^)11(^(^)3 < cellV X w|l(£2(fi))3, (40)

V х (a2(x,wi)wi) - V х {a2{x,W2)w2)\\(L2m3 ^ c7||™i - ^llw^fi)3* (41) где Ci > 0, сг — const, г = 5, 6. 7, a w, wi, w2 G ^¿(П)3.

Справедлива теорема. Теорема 2.2.2

Пусть w+ = tw2 + (1 — t)wi, t e (0,1) и wi,w2 б WjoW3, &сли для коэффициента а2 выполняются неравенства (40),(41) и справедливо представление а2{х, w+)w+ = ta2(x, w2)w2 + (1 - t)a2(x, w{)wi, то при / € L2(Q), F E (L2(Q,))3 функционал F2, определённый на И^о^)3, является непрерывным, коэрцитивным и строго выпуклым.

Обобщённая формулировка для задачи (37) получается с помощью дифференцирования функционала F2. При выполнении (40) и неравенства iiv х (v\a2)w,u)u)\\Lm < с8|ии,21о(п)з, где v'(a2. w, и) = {д(а2(x)w)wl) jди3}г^=1у2^ w,u G И^ДП)3, с8 = const, доказывается лемма 2.2.2 о том, что градиент функционала F2, определяемый равенством

A2w,u){L2{n)) з = J [(V'W-f)V-u+(Vx(a2(x,w)w)-F)-Vx(v'^ n задаёт отображение из в

При выполнении условий теоремы 2.2.2 и леммы 2.2.2 обобщённая формулировка для задачи (37) в виде уравнения

A2w,U)=Q, WeW^nf, имеет единственное решение w Е W^o^)3- Кроме того, при каждом п это уравнение для Vun £ (ЭД^о)".^)3 имеет единственное галёркинское приближение wn Е (И/21)0)п(^)3

Если в задаче (32) / = 0, то для её эффективного решения в полученной обобщённой формулировке рекомендуется использовать удовлетворяющий краевым условиям базис из пространства U(fl,).

В § 2.3 представлены проекционно-сеточные схемы для решения нелинейных задач вида (31),(32) относительно вектор-функций, использующие предложенные в п.2.1.3 функции формы из пространства У(П). Основой для построения этих схем является подход Бубнова-Галёркина.

Здесь строится проекционно-сеточная схема относительно вектор-функции для решения нелинейной краевой задачи с внутренней границей: V х ^ = 0, хеП] о, V х ¿г = Р'(х), у = г, хеп'-, (42) п-У] — 0, [пхг] = 0, хедП; V = г = 0, ж € г0.

Предполагается, что граница дО, является общей для областей П и ГУ, причём П' окружает П так, что Го является удалённой от частью границы 80,'. Важно отметить, что область может быть составлена из нескольких несвязных областей. Функция р задана либо как р = р(х,У), либо как р = р{х.Ё). Через [п-У] и [п х 7>] обозначена разность пределов соответствующих функций при стремлении точки х на границу дО, изнутри и извне области П. Относительно нелинейного коэффициента задачи предполагается, что выполняются обычные неравенства.

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

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

Третья глава диссертации посвящена методам и алгоритмам, существенно повышающим эффективность решения как прямых, так и обратных нелинейных задач магнитостатики. Среди публикаций на эту тему отметим работы Н.И.Дойникова, С.Е.Сычевского, Е.А.Ламзина с соавторами, J.С.Coulomb с соавторами, J.Simkin, C.W.Trowbridge с соавторами, M.J.Friedman, J.E.Pasciak, I.D.Mayergoyz с соавторами, Е.П.Жидкова, С.Б.Ворожцова, П.Г.Акишина, Б.Н.Хоромского, Э.А.Айряна с соавторами, А.Г.Дайковского с соавторами, S. Russenschuck с соавторами и других исследователей.

В § 3.1 предлагается новый метод вычисления функции на липшицевой границе трёхмерного тела по заданному градиенту. В качестве, примера применения метода приведём следующую нелинейную краевую задачу. Пусть область О' с R3 ограничена липшицевой границей и О' = n(Jf2i, где области Q и разделяются внутренней границей Г С О,', которая также является липшицевой. Требуется найти функции удовлетворяющие уравнениям

3 д cii(x, и, Vu)) + а(х, и, Vii) = 0, x G П; г=1 dXi

43) ij = 1 J i=1

3 3 Q u — ui = Q, ^(аг(х,и,¥и)- biUi) • ai = 0, xGT] г~ 1 j=1 3

3 а 3 и = й(х), х € Ш\Г; (JJ^jWq--= щ(х), х 6 <ЭПДГ, где Vu - градиент функции и, а о^о^ - компоненты векторов внешних нормалей к соответствующим областям. При этом равенства на границе понимаются в смысле следов соответствующих функций. Предполагается, что все коэффициенты в уравнениях вещественные и вместе с заданными функциями удовлетворяют условиям, необходимым для однозначной разрешимости задачи. Если в области f2i известно частное решение - достаточно гладкая вектор-функция G, компоненты которой G{ (i = 1, 2, 3) удовлетворяют уравнению то можно ввести функцию ф, такую, что \7ф = С в и она будет являться частным решением уравнения в области Г21. Тогда задачу (38) можно сформулировать относительно функций и в области Г2 и щ = щ — ф в области . Причём в области получается однородное уравнение для —. —* « функции щ. Поэтому если вектор (7 вычисляется более точно и'быстро, чем функция / в области и при этом упрощается процесс построения трёхмерной сетки, то переформулировка задачи для неизвестных и,й\ является оправданной. Однако точность решения новой задачи зависит от точности вычисления функции ф, которая входит в условия на границе Г.

Для определения функции ф введём пространство функций

Г) = {и:ие ь2(Т)ути е (Ь2(Г))2}, где в локальной правой декартовой системе координат с ортами 5, t, п, направленными вдоль касательной, бинормали и внешней нормали,

7тп = вди/дв + ¿Зи/дг.

В ]¥2(Г) выделяется класс функций И/21о(Г); для которых выполняется хотя бы одно из условий

J udS — 0 или J udl = 0.

Го la здесь Го - часть границы Г, a Iq - некоторая незамкнутая кривая на Г. Для функций u,v Е ^2",о (Г) можно ввести скалярное произведение и соответствующую норму

Ы, ^V2i0(r) = (VTU, Vr,)(L2(r))2, |М1ж,уг) = (0> u)wi0(T))1/2■

Эта норма эквивалентна норме пространства W^iT). Функция ф находится из обобщённого уравнения j Утф ■ VTvdS + Jp(G, \Чгф\){Утф - GT) ■ VTvdS = j GT ■ VTvdS,

Ti Г2 Г: где v E И^оСГ), Ti, Г2 - непересекающиеся составные части границы Г, определяемые в зависимости от поведения G Е К — {Чф : ф Е -^(Г), V</> Е (Ь2(Г))3)}. Здесь

Для исследования этого уравнения введём следующие обозначения. Зададим операторы В\, В2 с помощью равенств

BiUi,vi)wi^Tl) = J VTwi • VTvidS, m,vi E WlQ{Ti);

Г1

B2u2,v2)wio{r2) = Jp(G,\WTu2\)(VTu2 - GT) ■ VTv2dS, u2,v2 E W^Q^). r2

Очевидно, что Вг : WJqQTi) ^оО^) и В2 : Г2) И^0(Г2). Свойства оператора В2 сформулированы в следующей лемме. Лемма 3.1.1 Оператор В2 : W2 0(T2) —> W2 0(r2) является:

1) сильно монотонным, то есть выполняется неравенство

В2и ~ В2у, и - у)що{Г2] > \\и - и>у е ^2,о(г2);

2) липшиц-непрерывным, то есть выполняется неравенство

IIВ2и - В2у\\що(г2) < сг7\\и - и,уе И^2,о(г2), где 0-7 = 2 + тах(1/|£п|); бе К.

На основе леммы 3.1.1 доказывается, что оператор В : Н —» Н, где Н = ЖУГ) ГК^гО^) х ^(Гг)}, задаваемый равенством

Ви,у)н = (Вхиъи^щ^^ + (В2и 2,У2)\У10(Г2) для и = (и1,и2), v = (у\, у2) е н, является сильно монотонным и липшиц-непрерывным. Из общих теорем об уравнениях с такими операторами сразу же следует однозначная разрешимость обобщённого уравнения задачи и сходимость галёркинских приближений к точному решению. Кроме того, известен итерационный процесс, с помощью которого можно получить это решение. Предложенный метод вычисления функции на липшицевой границе трёхмерного тела по заданному градиенту по сравнению с "классическим"подходом [92] даёт приближения, сходящиеся более быстро к точному решению, поскольку использует двумерное интегрирование, а не интегрирование по отрезкам или кривым второго порядка. Кроме того, при решении системы уравнений не происходит накапливания ошибки, если градиент задан не точно.

В этом же параграфе приводится алгоритм вычисления функции на липшицевых границах несвязных областей по заданному градиенту для задач магнитостатики. Скалярный потенциал проводников с током (рБ предлагается определять на границах несвязных областей Г&, к = 1, 2,., т, из уравнений

Ус/^Я^Ог), хеТк, к = 1,2,.,т, (44) с условием ^(ук) = > 2//е ^ Гд;, (р3(уо) - известная величина, где г/о может не принадлежать границам Г^, А: = 1, 2,., т. В общем случае точки у к, к — 1,2связаны с точкой уо ломаными = 1,2, .,т, а потенциал на этих ломаных находится с помощью адаптивного алгоритма из обобщённых уравнений

Г д^ ду [ дхк \ ду

Ьк Ьк где у - достаточно гладкая функция, Хк ~ функция, обеспечивающая выполнение условия в точке уо, вектор ?{х) задаёт направление ломаной Ьк- Обобщённые линейные уравнения, соответствующие (44), будут иметь вид

J Ут</-Чтус13 = + А: = 1,2, .,ттг, (45)

Г/С Г*; где функции у, имеют тот же смысл, что и и, Хк в предыдущем уравнении соответственно.

Для решения уравнений (45) адаптивным методом конечных элементов рассматриваются четыре типа четырёхугольных элементов. При построении конечно-элементной системы используются лагранжевы элементы с 9-ю узлами (тип 2) и серендиповы элементы с 17-ю узлами (тип 3). Для вспомогательных вычислений применяются лагранжевы элементы с 6-ю (тип 1) и 15-ю (тип 4) узлами. Доказана лемма 3.1.3 об оценках функции интерполяпта через старшие производные известной вектор-функции. Эффективность адаптивного алгоритма демонстрируется на модельном примере и на вычислении скалярного потенциала обмотки для большого соленоидального магнита ЬЗ [8].

§ 3.2 диссертации посвящён алгоритмам с контролем точности для решения задач магнитостатики. Разработки алгоритмов контроля точности для проекционно-сеточных методов содержатся в многочисленных работах. Особенно актуальным этот вопрос является для численного решения трёхмерных эллиптических задач с сильно меняющимися решениями, в расчётной области которых содержатся подобласти с неодинаковыми дифференциальными уравнениями. При этом сами подобласти могут иметь сложную конфигурацию. Здесь можно отметить книги [16, 7, 111, 112, 21, 44, 22, 23, 25], в которых рассматриваются вопросы, возникающие при построении такого типа алгоритмов, и есть ссылки на оригинальные публикации. Различные алгоритмы контроля точности используются в некоторых коммерческих системах автоматизированного проектирования при решении линейных и нелинейных краевых задач с эллиптическими дифференциальными операторами второго порядка, например, [42, 44, 46, 45]. Все апостериорные методы анализа ошибки проекционно-сеточных методов можно разделить на три основные группы [109, 110, 39]: 1) методы, основанные на решении задачи в дополняющих друг друга формулировках, когда точное решение принимает промежуточное значение между приближёнными решениями одной и другой формулировки; 2) методы с многократным решением задачи; 3) методы с локальным оцениванием ошибки. Наиболее экономичной по вычислительным затратам считается третья группа методов, которая и развивается в последнее время особенно интенсивно. В этой группе выделяются методы с целенаправленной /гр-адаптивностью [61], а также метод суперсходящегося покусочного восстановления (superconvergence patch recovery method, [113]). Об актуальности вопросов, рассматриваемых в этой главе, свидетельствует практически постоянное представление докладов по контролю точности расчётов в задачах магнитостатики на регулярных конференциях Сот-pumag, посвящённых моделированию магнитных полей и смежным темам. В отличие от известных методов для этого класса задач, в третьей главе предлагается комплексный подход, объединяющий на основе единых характеристик как тестирование измеренных функций магнитного поля спектрометров для физических экспериментов, так и тестирование магнитного поля, которое получается в результате расчётов.

В § 3.2 приводятся формулы, характеризующие численную погрешность в зависимости от невязки уравнений для операторов дивергенции и ротора. Пусть область О С К3 разбита стандартным образом на некоторые конечные элементы гиг-, г = 1,2,., ТУ, и на каждом г-ом элементе определены функции формы Л^Дж), Л^Дж/с) = 1; хк £ Щ- Через и)к обозначим носитель к-той базисной функции к — и зиРРМг,к,

Хк&Вг а через (Со°(щ))3 - множество финитных бесконечно-дифференцируемых вектор-функций с носителями на щ. Справедлива следующая лемма. Лемма 3.2.1

Пусть для V Е (Со°(а>а;))3 выполняются неравенства

IV • ф) - V • у(у)I < сА IV х - V х у {у)I < с2сг0, г(?е х,у Е сЬк, ¿о = тах |ж — сх,с2 - положительные постоянные, не зависящие от х и у. Тогда при х Е ¡2>к выполняется равенство у{х) = (V • у(х))То(х) + (V х у(х)) х /о(ж) + Л0V, здесь

Нх) = — / ^у, |/о(ж)| < ¿0, < к 3 при этом с3 = тах{с1,с2}; г = ^ гк(хк - ук), г =\х -у\. к=1

Справедлива также лемма о представлении вектор-функции V Е (С°°(и1/;;))3 на элементе и^. Лемма 3.2.2

Пусть для V Е (С°°(шк))3 выполняются неравенства

У • у(х) - V • у{у)| < IV х у(х) - V х гГ(т/)| < с5(1, где х:у Е (¿к, ^ ~ шах \х — х\, С4, £5 - положительные постоянные, не зависящие от х и у. Тогда при х Е и!к выполняется равенство и{х) = (V ■ у(х))1(х) + (V х гГ(ж)) х 1{х) + Ру + Ш, здесь

Шк 3

РУ\ < тах(У^ Ы)Ы> |Яг;| < М2, с6 = тах{с4, с5}. гС— 1

Леммы 3.2.1 и 3.2.2 можно использовать для оценки погрешности приближённого решения краевых задач с операторами дивергенция и ротор.

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

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

77 (1®В т) = (| J V • + \ ! V х 1^Вс11и\)/{\<ш\В0), ги и> где и) - некоторый конечный элемент, образованный из узлов сетки, на которых известны данные измерений, |ги| - объем ги, а Во -величина основной компоненты поля в центре магнита. Через обозначен оператор восполнения. В некоторых случаях он может быть равен оператору 1е \ задаваемому формулой mi

B^B.Jvfw, г=1 где индекс I - степень интерполяционного многочлена, то/ - число используемых при интерполяции узлов сетки, Вг - значение вектора В в г-м узле, N^l\x) - г-я базисная функция интерполяционного многочлена.

С помощью характеристики г], во-первых, можно определить ошибки в измерениях поля, а во-вторых, можно установить, на какой сетке необходимо проводить измерения, чтобы ошибка функции магнитного поля была минимальной. В качестве критерия следует использовать неравенство r](I^B,w) <С -hl+1 -\Dl+1B{x)\, xew, где С = const, С > 0, h - диаметр элемента w, \m\ = mi -Ь Ш2 + тоз.

В § 3.3 представлены методы решения задачи формирования однородного магнитного поля, которую можно отнести к классу обратных задач. Общий подход к их решению разработан А.Н.Тихоновым и его школой. Здесь следует отметить работы П.Н.Вабищевича [114], а для задач магнитостатики - работы Е.А.Ламзина, С.Е.Сычевского с соавторами [115, 116]. В этом разделе сформулирован итерационный процесс уточнения части поверхности ферромагнетика для получения высокооднородного магнитного поля (теорема 3.3.1) и даны примеры решения задач с его помощью.

В § 3.4 приводится алгоритм построения компьютерных моделей анализирующих магнитов спектрометров с контролем точности вычислений. Главными требованиями к таким магнитам обычно являются: величина поворотной силы, объём рабочей области, эксплуатационные условия, материал проводника и тип стали. Процесс построения компьютерной модели дипольного магнита, основанный на этих входных данных, можно описать с помощью следующего алгоритма:

1) оценка ампер-витков, длины магнита и толщины ярма по аналитическим формулам;

2) построение компьютерной модели обмотки и итерационный процесс формирования геометрии ненасыщенного магнитного ярма;

3) анализ точности разработанной компьютерной модели, который включает в себя проверку правильности вычисления поля от обмотки и анализ адекватности расчётной сетки;

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

Четвертая глава диссертации посвящена результатам применения методов, описанных в предыдущей главе, к построению компьютерных моделей и расчётов полей магнитных систем для установки ALICE [8]-[11]. (CERN, г. Женева) и проекта установки для эксперимента с поляризованной мишенью в ИТЭФ [12] (г. Москва).

В § 4.1 приводятся результаты построения компьютерных моделей для большого соленоидального магнита L3, проекта сверхпроводящего дипольного магнита и первого проекта "теплого"дипольного магнита для эксперимента ALICE. В частности, для модели L3 представлены результаты сравнения расчётов, выполненных по программе MSFE3D [167], с картой поля магнита, а также с результатами расчётов по другим программам.

В § 4.2 представлены результаты построения компьютерных моделей дипольных магнитов серии VULCAN с конусной конфигурацией ферромагнетика для эксперимента ALICE. Здесь приводятся разработанные модели обмоток, укладываемых на боковую поверхность прямого усечённого конуса, результаты использования различных характеристик для оценки точности вычисления поля в апертурах магнитов, а также другие расчёты.

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

К диссертационной работе имеются три приложения.

Приложение 1 содержит анализ функций магнитного поля, основанных на данных измерений для эксперимента ЭКСЧАРМ [120]. Рассматриваются 4 различные функции магнитного поля, интерполирующие данные с сетки измерений с помощью произведения интерполяционных многочленов Лагранжа первой, второй третьей и четвертой степеней соответственно по каждой пространственной переменной. Для их сравнения используется характеристика t](I^B,wq), 1 < к < 4, введённая в § 3.2.

В Приложении 2 приводятся результаты расчётов магнитного поля для модели соленоида проекта модифицированного бетатрона [13], разработанного в Лаборатории Ядерных Проблем ОИЯИ.

В Приложении 3 представлены результаты расчётов магнитного поля для модели поворотного магнита М-90-05 ускорительного комплекса У400-У400М, предназначенного для проекта DRIBs [14].

Основные результаты диссертационной работы докладывались на международных совещаниях по проблемам математического моделирования, программированию и математическим методам решения физических задач (Дубна, 1983, 1993), на X, XI, XV Всесоюзных и Всероссийском совещании по ускорителям заряженных частиц (Дубна,

1986, 1988, Протвино, 1996), на II Республиканской конференции "Интегральные уравнения в прикладном моделировании"(Киев, 1986), на Всесоюзных конференциях по вычислительной физике и математическому моделированию (Волгоград, 1988, 1989), на Рабочем совещании "Методы и программы расчёта магнитных полей" (Протвино, 1989), на международных конференциях по дифференциальным уравнениям EQUADIFF -7, 8 (Прага, 1989, Братислава, 1993), на Международной конференции по численному анализу ISNA'92 (Прага, 1992), на конференции по ускорителям заряженных частиц РАС-95 (Dallas, 1995), на рабочих совещаниях коллабораций ALICE (CERN, 1996), EXCHARM (Дубна, 1995, 1997) и PANDA (FZ-Juelich, 2004), на конференции-семинаре "Математические модели сложных систем" (Тверь, ТвГУ, 1999), на всероссийских конференциях по проблемам физики, .химии, математики, информатики и методики преподавания (Москва, РУДН, 1999, 2000, 2002, 2006), на международных конференциях "Тихонов и современная математика" (Москва, МГУ, 2006), "Математическое моделирование и вычислительная физика" (Дубна, 2009), на семинаре кафедры Вычислительных методов (ВМК, МГУ, 2009), на семинаре Научно-исследовательского вычислительного отдела НТЦ «Синтез» (НИИЭФА, С.-Петербург, 2009), а также на Ученых советах ОИЯИ и семинарах по вычислительной и прикладной математике, вычислительной физике ЛВТА, ЛИТ ОИЯИ.

Основные результаты, вошедшие в диссертацию, опубликованы в 28 печатных работах, в том числе в 9 публикациях в журналах, рекомендованных Высшей аттестационной комиссией ("Дифференциальные уравнения-[85], "Математическое моделирование-[70], "Журнал вычислительной математики и математической физики-[69], "Вестник Российского Университета Дружбы Народов-[121, 78, 122, 123, 124], "Nuclear Instruments and Methods-[13]), в журнале "Краткие сообщения ОИЯИ-[125], а также в 18 публикациях в трудах Всесоюзной [127] и Международных конференций [86, 132, 133, 79, 145], в сообщениях ОИЯИ [71, 134, 135, 136, 137, 138, 139], а также в других научных изданиях [140, 141, 142, 80, 143, 144]. Некоторые из работ автора, на основе которых написана диссертация, представлены в электронной библиотеке публикаций ОИЯИ [146], в электронных библиотеках институтов КЕК (Япония) и CERN (г.Женева) [147].

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

ЗАКЛЮЧЕНИЕ

Приведём основные результаты диссертации.

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

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

3. Разработаны методы и алгоритмы решения трёхмерных нелинейных задач магнитостатики:

- проекционно-сеточный метод и адаптивный алгоритм вычисления скалярного потенциала, создаваемого проводниками с постоянным током, на границах несвязных областей вне проводников;

- алгоритм построения единых характеристик для контроля точности вычисленных и измеренных магнитных полей ;

- итерационный метод формирования высокооднородного магнитного поля за счёт изменения объема, занимаемого ферромагнетиком;

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

4. На основе созданных алгоритмов и программного обеспечения решён ряд трёхмерных задач с учётом нелинейного влияния ферромагнитных частей магнитных систем. В некоторых задачах использовались характеристики точности вычислений. Разработаны компьютерные модели и проведены расчёты:

- большого соленоидального магнита L3 установки ALICE на LHC и проектов ОИЯИ дипольных магнитов для этой установки: сверхпроводящего и "тёплых"вариантов, в том числе с конусной конфигурацией магнитопровода;

- спектрометрического магнита с поляризующими наконечниками для проекта эксперимента с поляризованной мишенью на ускорителе в ИТЭФ;

- соленоида для проекта модифицированного бетатрона ОИЯИ.

Автор выражает глубокую благодарность своим научным и консультантам: заслуженному деятелю науки РФ Жидкову Е.П. доктору физико-математических наук Иванову В.В. за поддержку и помощь на всех этапах работы над диссертацией. Автор благодарен заслуженному деятелю науки РФ Пузынину И.В., профессору Позе Р.Г., другим членам дирекции ЛИТ (ЛВТА) за предоставление благоприятных условий работы над диссертацией.

Автор благодарен доктору физико-математических наук Акишину П.Г. за плодотворное сотрудничество и полезное обсуждение результатов диссертации. Автор глубоко благодарен кандидату физико-математических наук Юлдашевой М.Б. за помощь и поддержку, которую трудно переоценить.

Автор благодарен докторам физико-математических наук Водопьянову A.C., Иванченко И.М., Маторе И.М., Рыльцову В.В.(ИТЭФ), за плодотворное сотрудничество и поддержку, кандидатам физико-математических наук Айряну Э.А., Федорову A.B. , Новицкому В.В., Сидорину А.О., Шишову Ю.А., Юдину И.П., за сотрудничество и полезные обсуждения, кандидатам физико-математических наук Н.Н.Карпенко, Т.А.Стриж, Ж.Ж.Мусульманбекову, докторам М. Pavlus (University of Presov, Slovakia), A. Polansky, профессорам M.Slodicka (Ghent university, Belgium), M.Gregus (University of Bratislava, Slovakia), за полезные обсуждения и поддержку. Автор благодарен коллегам из коллабораций ЭКСЧАРМ, ALICE, профессору J. Ritman (FZ-Juelich, Germany) за сотрудничество, сотрудникам Научного Отдела Вычислительной Физики ЛИТ за полезные обсуждения и поддержку, а также многим сотрудникам ЛИТ (JIBTA) за помощь во время работы над диссертацией. Автор также выражает особую благодарность Г.Г.Гульбекяну и В.В.Башевому за сотрудничество.

1. Крылов А.Н. Лекции о приближенных вычислениях. 5-е издание, М.-Л., ГИТТЛ, 1950, 400с.

2. Канторович Л.В., Крылов В.И. Приближенные методы высшего анализа. М.-Л., ГИТТЛ, 1952, 695с.

3. Березин И.С., Жидков Н.П. Методы вычислений, т. 1,2, М., Физматгиз, 1966.

4. Марчук Г.И. Методы вычислительной математики. М. Наука, 1989, 608с.

5. Самарский A.A. Теория разностных схем. М., Наука, 1983.

6. Самарский A.A., Гулин A.B. Устойчивость разностных схем. М., Наука, 1973.

7. Бахвалов Н.С. Численные методы. М., Наука, т.1, 1975.

8. ALICE collaboration.(Ahmad N., . Vodopianov A.S., .YuldashevO.I., Yuldahseva M.B. et al.) ALICE techinal proposal. CERN-LHCC-95-71, 1995, 237p.

9. ALICE collaboration (Beole S., . Vodopianov A.S., . Iouldachev O.I., . Yuldahseva M.B. et al.) The Forward Muon Spectrometer of ALICE (Addendum to the ALICE technical proposal), CERN/LHCC 96-32,1996, 77p.

10. ALICE collaboration. (Beole S., .Vodopianov A.S., . Iouldachev O.I., . Yuldahseva M.B. et al.) ALICE technical design report: detector for high momentum PID. CERN-LHCC-98-19, 1998, 198p.

11. ALICE collaboration. (Dellacasa G., . Vodopianov A.S., .Iouldachev O.I., . Yuldahseva M.B. et al.) ALICE technical design report: photon multiplicity detector (PMD). CERN-LHCC-99-32, 1999, 134p.

12. Алексеев И.Г., . Рыльцов В.В. и др. Исследование поляризованных параметров в бинарных реакциях рождения странных частиц ттр —> КА(Т<) на ускорителе ИТЭФ (Предложение эксперимента), ГНЦ РФ ИТЭФ, 41-97, Москва, 1997.

13. Ivanov A., Ivashkevich S., Korotaev Yu., Meshkov I.,., Yuldasheva M., Yuldashev 0. Focusing system of the modified betatron: design, technology, manufacturing and test. NIM, A, v.441, 2000, p.262-266.

14. Gulbekian G.G., Oganessian Y.T. "DRIBs": The Dubna Project for Ra-diactive Ion Beams. Nuclear Shell 50 Years: Int. Conf. on Nucl. Phys. (Ed. by Y.T.Oganessian), Singapore etc., World Sci., 2000.

15. Тихонов A.H., Арсенин В.А. Методы решения некорректных задач. М., Наука, 1979, с.18

16. Михлин С.Г., Смолицкий X.JI. Приближенные методы решения дифференциальных и интегральных уравнений. М., Наука, 1965.

17. Красносельский М.А., Вайнокко Г.М., Забрейко П.П. и др. Приближенное решение операторных уравнений. М., Наука, 1969.

18. Сьярле Ф. Метод конечных элементов для эллиптических задач. М., Мир, 1980, 512с.

19. Оганесян JI.A., Руховец JI.A. Вариационно-разностные методы решения эллиптических уравнений. Ереван, Изд-во АН Арм.ССР, 1979.22