Содержание к диссертации
Введение
Глава 1. Методы решения задач магнитостатики с использованием МКЭ и МГЭ 10
1.1. Математические модели магнитостатики 10
1.1.1. Модель с полным скалярным потенциалом 11
1.1.2. Модель с неполным скалярным потенциалом 12
1.1.3. Модель с двумя потенциалами 13
1.1.4. Обобщенная модель 14
1.2. Метод граничных элементов 15
1.2.1. Вариационные постановки для краевых задач 17
1.2.2. Дискретизация вариационных уравнений 20
1.2.3. Учет симметрии в МГЭ 23
1.2.4. Вычисление решения в МГЭ 23
1.3. Метод конечных элементов 24
1.3.1. Вариационная постановка и дискретизация 24
1.3.2. Вариационная постановка для нескольких подобластей
1.4. Вычисление производных решения по параметрам геометрии 26
1.5. Метод декомпозиции на подобласти
1.5.1. Описание FETI 29
1.5.2. Предобусловливание 33
Глава 2. Вычисление потенциалов 35
2.1. Методы вычисления несобственных интегралов 35
2.2. Аналитическое вычисление потенциалов
2.2.1. Вычисление потенциала простого слоя 37
2.2.2. Вычисление потенциала двойного слоя 40
2.2.3. Вычисление потенциала простого слоя от вектора 41
2.2.4. Вычисление потенциала двойного слоя от вектора 42
2.3. Учет токовых обмоток 44
2.3.1. Вычисление напряженности поля токовых обмоток 44
2.3.2. Вычисление вектор-потенциала поля токовых обмоток з
2.4. Использование мультипольного разложения 48
2.5. Выводы 52
Глава 3. Описание программного комплекса Quasar 53
3.1. Структура программного комплекса Quasar 53
3.2. Автоматизация построения согласованных базисов высокого порядка в методе конечных элементов 56
3.3. Вычисление локальных матриц в комплексе Quasar 60
3.4. Выводы 61
Глава 4. Верификация 62
4.1. Сравнение результатов моделирования с измерениями 62
4.2. Сравнение с комплексом TELMA при моделировании H-образной балки 65
4.3. Сравнение реализованных методов между собой, на примере T-образной балки 67
4.4. Оценка эффективности FMM на задаче с большим количеством объектов 71
4.5. Проверка процесса подбора геометрии на тесте с известным ответом 72
4.6. Выводы 76
Глава 5. Решение практических задач 77
5.1. Моделирование возмущений магнитного поля внутри помещений 77
5.2. Оптимизация геометрии ускорительных магнитов 81
5.3. Выводы 90
Заключение 91
Список литературы
- Вариационные постановки для краевых задач
- Вычисление потенциала двойного слоя
- Автоматизация построения согласованных базисов высокого порядка в методе конечных элементов
- Оптимизация геометрии ускорительных магнитов
Введение к работе
Актуальность темы исследования следует из необходимости высокоэффективных методов численного моделирования при проектировании технических устройств. Сростом сложности технических устройств растути требования к точности требуемых для их проектирования расчетов. Так, при проектировании современных ускорительных магнитов зачастую требуется вычисление магнитного поля с точностью до сотых долей процента. Для получения требуемой точности с использованием метода конечных элементов требуются подробные сетки очень высокого качества. Особенно остро эта проблема проявляется в задаче автоматической оптимизации формы ускорительных магнитов, когда требуется постоянно перестраивать сетку и вычислять производные минимизируемого функционала в зависимости от оптимизируемых параметров геометрии магнита. Сохранять в этих условиях качество конечноэлементной сетки, с учетом того, что сильные изменения сетки будут крайне негативно отражаться на процессе минимизации, становится весьма нетривиальной задачей.
Применение подхода с совместным использованием методов конечных и граничных элементов позволяет избавиться от необходимости построения сетки в большой части расчетной области, что дает возможность значительно упростить построение сетки. Помимо этого, применение метода граничных элементов позволяет получить значительно более гладкое решение, менее чувствительное к особенностям сетки, что улучшает сходимость методов оптимизации.
Используемые в настоящее время для расчетов ускорительных магнитов программные комплексы, такие как MERMAID, Opera3D, MASTAC и ANSYS, не используют метод граничных элементов. В разработанном CERN программном комплексе ROXIE применяется совместный метод конечных и граничных элементов, однако используется менее эффективная постановка с векторным магнитным потенциалом.
Цель работы состоит в разработке и реализации методов численного моделирования для эффективного решения задач магнитостатики, возникающих при проектировании технических устройств. В соответствии с поставленной целью предусмотрено решение следующих задач.
-
Разработать вычислительную схему с возможностью декомпозиции расчетной области на подобласти, позволяющую выбирать способ аппроксимации в каждой подобласти, с возможностью использования коэффициента магнитной проницаемости, зависящего от магнитного поля.
-
Разработать эффективные методы вычисления магнитного поля токовых обмоток в однородном пространстве.
-
Исследовать возможность моделирования возмущений магнитного поля внутри помещений.
4. Исследовать эффективность различных вариантов решения задачи оптимизации геометрии ускорительных магнитов
Научная новизна работы:
-
Разработана вычислительная схема для совместного использования конечных и граничных элементов при численном решении задач магнитостатики на основе модели с полным и неполным скалярными потенциалами.
-
Разработаны и реализованы эффективные вычислительные схемы метода граничных элементов на основе мультипольного разложения.
-
Разработан исследовательский программный комплекс Quasar, в котором реализована технология автоматического согласования базисных функций и расчета локальных матриц, для которой необходимо задать только выражения для локальных базисных функций.
Методы исследования основаны на использовании метода конечных и граничных элементов для решения трехмерных нелинейных задач магнитостатики, быстрого мультипольного метода для ускорения метода граничных элементов, методов многомерного нелинейного программирования для оптимизации геометрии ускорительных магнитов.
На защиту выносятся:
-
Вычислительная схема для решения трехмерных нелинейных задач магнитостатики, допускающая совместное использование конечных элементов в подобластях с полным скалярным магнитным потенциалом и граничных элементов в подобластях с неполным скалярным магнитным потенциалом, позволяющая учитывать нелинейные свойства ферромагнетиков и условия симметрии моделируемого процесса.
-
Результаты исследований эффективности мультипольной реализации метода граничных элементов, на примере решения задачи с большим количеством ферромагнитных объектов.
-
Эффективные алгоритмы вычисления векторного магнитного потенциала и индукции магнитного поля токовых обмоток в однородном пространстве.
-
Объектно ориентированный программный комплекс Quasar, реализующий все описанные в работе вычислительные схемы и алгоритмы.
Достоверность результатов обеспечивается корректным применением доказанных теоретических положений, апробированных методов, тестированием разработанных программ, верификацией результатов на модельных задачах, путем независимого решения их различными методами, сравнением решений задач с решениями других авторов, сравнением результатов моделирования с результатами экспериментальных электромагнитных измерений.
Практическая значимость и реализация результатов
Практическая значимость работы состоит в возможности использования разработанных математических моделей, численных методов и реализующего их комплекса программ Quasar при решении задач проектирования и автоматической оптимизации геометрии магнитов в ускорителях заряженных частиц.
Результаты работы были использованы при выполнении следующих работ:
Хоздоговорная работа «Конечноэлементные исследования магнитных полей косинусных магнитов» (2009 г., НИУ ИЯФ СО РАН).
Государственного контракта № 14.740.11.0709 (разработка комплекса программ для автоматизации 3D проектирования дипольных магнитов) от 12 октября 2010 г.
Договора с компанией Samsung по теме «Sophisticated Modeling of Indoor Magnetic Field Disturbance» (Сложное моделирование распределения геомагнитного поля в здании) за 2011-2012 гг.
Теоретическая значимость работы состоит в том, что исследована эффективность различных вычислительных схем с совместным использованием конечных и граничных элементов.
Личный вклад автора заключается в формулировке математической модели, разработкеиреализации вычислительных схем, получении аналитических выражений для используемых потенциалов, проведении вычислительных экспериментов, разработке комплекса программ Quasar для решения трехмерных нелинейных задач магнитостатики.
Апробация работы. Основные положения и результаты работы докладывались на Всероссийской научной конференции молодых ученых «Наука Технологии Инновации» (Новосибирск, 2006, 2010, 2012 и 2013), Всероссийской конференции по вычислительной математике (Новосибирск, 2009), XVII международной конференции «Математика. Компьютер. Образование» (Дубна, 2010), VII Всероссийской межвузовской конференции молодых ученых (Санкт-Петербург, 2010), Российской научно-технической конференции «Обработка информационных сигналов и математическое моделирование» (Новосибирск, 2012), I Всероссийском конгрессе молодых ученых (Санкт-Петербург, 2012), международной конференции «Актуальные проблемы электронного приборостроения» (Новосибирск, 2012 и 2014), международной конференции «Актуальные проблемы вычислительной и прикладной математики» (Новосибирск, 2014).
Публикации: основные положения диссертационной работы опубликованы в 12 работах, в том числе 4 статьи в изданиях, рекомендованных ВАК, 8 статей в сборниках трудов.
Структура и объем работы. Диссертация состоит из введения, пяти глав, заключения, списка использованной литературы, списка иллюстративного материала и приложений. Общий объем основной части составляет 103 страницы и включает 42 рисунка и 3 таблицы.
Вариационные постановки для краевых задач
Поскольку модель с неполным потенциалом обладает рядом серьезных недостатков, а использование полного возможно только в достаточно редко встречающихся задачах, была предложена комбинированная модель [18; 28; 30; 31; 80; 81], в которой используются сразу два потенциала. Полный потенциал обычно используется внутри металлов, а неполный в воздушной среде. Отметим, что для выполнения закона Ампера требуется, чтобы нельзя было обойти ненулевой ток по подобласти полного потенциала.
Таким образом, в областях с полным потенциалом должно выполняться уравнение (1.6), в областях с неполным потенциалом уравнение (1.10), на границах между областями с полным потенциалом условия (1.7), а на границах между областями с неполным потенциалом условия (1.11).
Осталось определить условия на границе между полным и неполным потенциалами. Обозначив эту границу через TR, условия связи можно записать в виде
Из условия (1.13) видно, что при вычислении Hs через закон Био-Савара-Лапласа (1.4) определить потенциал всюду непрерывным нельзя1. Введем функцию разрыва в так, чтобы на границе Гуд выполнялось равенство ит — UR = О. Тогда для выполнения условия (1.13) требуется, чтобы функция в удовлетворяла
Отметим, что впринципе возможно использовать модель с полностью непрерывным потенциалом [66],но для этого требуется выделить внешнее поле H s таким образом, чтобы на границе TR оно было равно нулю. следующему уравнению:
Пусть расчетная область П разбита на непересекающиеся подобласти Пг. Обозначим через Г границу Г , через Г границу между Г и ,а через щ, [м и Mj - соответственно потенциал, магнитную проницаемость и намагниченность в Пг. Для записи условий на границах областей введем операторы следа Дирихле 7о и следа Неймана 7І: У0и (х) = lim и (г), 7о : Н1 (ПА Я1/2 (ГЛ , 7Їи (х) = lim пг (х) gradw (г), 7! : я1 №) я"1/2 (г ), ГЄПІ;Г-4-ХЄГІ где пг - внешняя нормаль к области Пг, а Я1 (П{), Я1/2 (Гг) и Я"1/2 (Гг) - пространства Соболева [1; 82] с соответствующим нормами:
Будем полагать, что все М; являются постоянными и не равны нулю только во внутренних подобластях (этого достаточно для учета постоянных магнитов, встречающихся в задачах магнитостатики). Также будем полагать, что в подобластях с неполным потенциалом / являются константами. С учетом этих предположений, уравнения (1.6) и (1.10) принимают вид -div/i,grad (x) = 0, х є Пг с Пт, (1.16) -Аи{(х) = 0, ХЄПІСПД, (1.17) где функции щ Є Нг{Пг). Будем считать, что на внешнюю границу Q выходят только области с неполным потенциалом. Тогда условия на границах можно записать в виде
Первый интеграл в формуле (1.25) называется потенциалом простого слоя, а второй – потенциалом двойного слоя. Отметим, что с учетом (1.22) формула (1.25) остается верной и для случая неограниченной области.
При помощи формулы (1.25) можно вычислить решение в любой точке внутри области, если известны значения функции и потока на границе. Также из нее можно получить граничные интегральные уравнения, на которых основывается метод граничных элементов. 71 intM = a-f1 intu + K 7intw + D70V (1.30) где K - сопряженный оператор двойного слоя, а D - гиперсингулярный граничный интегральный оператор, который является ограниченным самосопряженным и полу-эллиптичным в пространстве Н1/2 (Г). Слагаемое а и появляется из-за того, что производная потенциала простого слоя имеет разрыв на границе.5 Преобразуем уравнение (1.30) к виду, удобному для нахождения следа Дирихле D u = ( I -K ) 71 inV (1.31) Рассмотрим вопрос существования решений уравнения (1.31). Для этого, подставив функцию щ = 1 в уравнения (1.28) и (1.31), получим следующие соотношения Kм0 = 0 Dw0 = 0. Можно показать, что ядро оператора D состоит только из постоянных функций. Следовательно, решение определено с точностью до произвольной константы, а для разрешимости уравнения (1.31) необходимо и достаточно, чтобы правая часть была ортогональна функции щ. С учетом сопряженности операторов K и K, получаем Заметим, что в случае внешней задачи оператор K меняет знак, следовательно, правая часть уравнения (1.31) всегда ортогональна функции щ. Соответственно, внешняя задача разрешима всегда. Причем, несмотря на то, что функция щ является решением уравнения (1.31) и в случае неограниченной области, она не удовлетворяет условию (1.22). Следовательно, решение внешней задачи единственно.
Для получения слабой формы уравнения (1.31) воспользуемся методом Га-леркина. Чтобы обеспечить единственность решения, можно искать его в подпространстве Я1/2 (Г), ортогональном щ. Однако, как показано в [82], вместо этого можно модифицировать билинейную форму оператора D так, чтобы она стала эллиптичной, и сделать вариационную постановку для пространства Я1/2 (Г)
Вычисление потенциала двойного слоя
Чтобы при использовании обычных методов численного интегрирования для вычисления сингулярных интегралов добиться достаточного уровня точности, требуется адаптивно сгущать сетку интегрирования к точке сингулярности. Аналогичное сгущение требуется даже если интеграл является собственным, но близким к сингулярному Это приводит к неприемлемому росту затрат времени на вычисление интегралов. Поэтому для вычисления таких интегралов используют специальные методы. Их можно разделить на две основных группы.
К первой группе относятся методы компенсации особенности с помощью специальной замены переменных. Такой подход был предложен в работе Duffy [39] и с тех пор активно развивается [32; 41; 49; 50; 79].
Второй подход заключается в выделении сингулярной особенности в отдельный, более простой, интеграл [47]. Подынтегральное выражение разбивается на две части, одна из которых является ограниченной и поэтому может быть проинтегрирована численно, а вторая не зависит от вида базисных функций и может быть проинтегрирована аналитически.
Рассмотрим вычисление потенциала простого слоя где if - базисная функция, являющаяся полиномом на элементе Г.
В случае, когда точка y принадлежит элементу Г, подынтегральная функция в этой точке является неограниченной и интеграл становится несобственным. Воспользуемся идеей, предложенной в [26], и выделим особенность в отдельный интеграл. Для этого добавим и вычтем из числителя функцию if в точке y. Разбив результат на два слагаемых, получаем а значит, подынтегральная функция в первом слагаемом будет ограниченной. Это, в свою очередь, позволяет применить для вычисления этого интеграла численное интегрирование. Второй интеграл не зависит от вида функции if и может быть вычислен аналитически.
Для вычисления потенциала двойного слоя можно применить аналогичный подход. Однако, из-за более высокого порядка особенности, в этом случае потребуется два слагаемых из ряда Тейлора. В результате получим
Ограниченность подынтегрального выражения в первом интеграле следует из соотношения ф (x) - ф (y) - (x - y) grad V (y) = О (\x - y2 ) . (2.6) Остальные интегралы, как и в случае потенциала простого слоя, могут быть вычислены аналитически.
В случае, когда точка не принадлежит к элементу Г, но находится достаточно близко (например, на соседнем элементе), интеграл становится близким к несобственному и численное интегрирование потенциалов напрямую дает большую погрешность. Для ее уменьшения можно также использовать предложенный выше подход. Однако, базисные функции обычно определены только на элементе Г, что не позволяет напрямую вычислять их в точке y. Доопределим их следующим образом (р(y) = (р (p) + (y - p) grad (р (p), (2.7) ф(y) = ф (p) + (y - p) grad V (p), (2.8) где p - проекция y на элемент Г. Легко проверить, что доопределенные таким образом базисные функции остаются гладкими, что позволяет использовать предложенный выше подход.
Будем считать, что область Г является плоской, тогда можно определить цилиндрическую систему координат (рисунок 2.1) следующим образом. За ноль возьмем точку y. Ось z направим перпендикулярно плоскости, в которой лежит область Г. Точка x в этой системе координат принимает вид (г, (p,z). Расстояние между точками x и y в можно записать следующим образом x - y = л/г2 + z2. (2.10)
Для того, чтобы в формуле (2.9) перейти к интегралу по границе области Г, представим подынтегральное выражение в виде ротора некоторой функции. С учетом того, что нормаль к области сонаправлена с осью z, а подынтегральное выражение не зависит от (р, получаем следующее уравнение
Далее будем предполагать, что область Г является многоугольником, и, следовательно, ее граница может быть представлена как набор отрезков. Рассмотрим вычисление интеграла (2.14) по отрезку. Повернем систему координат таким образом, чтобы при р = 0 касательная к прямой, на которой лежит отрезок, совпадала с вектором е . Тогда уравнение прямой можно записать следующим образом: LW
Чтобы привести полученное выражение к более удобному виду, определим еще одну систему координат (U,V,W).За ноль возьмем проекцию точки у на плоскость, в которой лежит область Г. Ось w направим перпендикулярно этой плоскости, ось и - перпендикулярно рассматриваемому отрезку, а ось v - параллельно ему (см. рисунок 2.2). В этой системе координат г о = щ, tg if = v/щ, h - расстояние от точки y начала координат Таким образом, в новой системе координат выражение (2.16) с учетом формулы (2.17) принимает следующий вид
Автоматизация построения согласованных базисов высокого порядка в методе конечных элементов
На каждой сетке для решения задачи использовалось три базиса: линейный, квадратичный и кубический. Также, на двух первых сетках было получено решение с использованием граничных элементов в воздушной области (при этом она считалась неограниченной). За точное было принято решение с использованием Для тестирования постановок с использованием различных потенциалов были проведены расчеты методом конечных элементов с использованием только полного потенциала, только неполного потенциала и двух потенциалов, полного внутри балки и неполного в воздушной среде. На данной задаче решения, полученные с использованием этих моделей, совпали, с точностью до погрешности решения СЛАУ. Поэтому все результаты будем приводить только для постановки с двумя потенциалами.
На рисунке 4.12 изображено поведение относительной погрешности с использованием различных методов на грубой сетке. При этом, для вычисления поля H в методе конечных элементов использовались согласованные результанты [18]. Как видно из рисунка, метод конечных элементов при использовании линейного базиса дает погрешность порядка 10%. В то же время, с использованием граничных элементов в воздушной среде, погрешность на линейном базисе не превышает
Относительная погрешность на разных сетках, без использования согласованных результантов На рисунке 4.13 показано поведение погрешности на различных сетках. Отличие порядков сходимости от теоретических может быть обусловлено использованием согласованных результантов для сглаживания решения в методе конечных элементов и недостаточной точностью эталонного решения для метода граничных элементов. На рисунке 4.14 изображено поведение погрешности, если не использовать согласованные результанты. Из рисунка видно, что, хотя поведение погрешности и становится ближе к теоретически ожидаемому, сама погрешность решения становится больше, особенно при использовании линейных элементов.
Для тестирования эффективности предложенной схемы решим задачу о вычислении искажений магнитного поля Земли металлическими объектами в стенах помещения. Расположение металлических объектов изображено на рисунке 4.15
В модели задано 319 металлический прутьев с квадратным сечением 22см, относительная магнитная проницаемость металла задана равной 1000. При решении этой задачи даже на очень грубой сетке обычным МГЭ затраты времени и памяти становятся значительно больше, чем при использовании ускоренного подхода. Так, прямой метод требует только на сборку матриц более часа времени и более 4ГБ памяти. ИспользованиеМКЭдля решения этой задачи также затруднительно – поскольку требует построения сетки в воздухе между объектами. Применение же быстрого МГЭ позволяет решить эту задачу даже на относительно подробной сетке за приемлемое время. Как видно из таблицы 4.1, для быстрого МГЭ рост затрат времен и памяти при дроблении сетки близок к линейному, тогда как у обычного МГЭ он является квадратичным.
Для тестирования подбора геометрии будем использовать задачу определения параметров металлического прямоугольного параллелепипеда в постоянном магнитном поле. Параллелепипед задается при помощи 6 параметров. Первые три параметра - х, у и z, задают центр параллелепипеда (все координаты от —0.9 до 0.9), а параметры sx, sy и sz - отвечают за расстояние между центром и гранями (от 0.1 до 1.0 по каждому направлению). Коэффициент магнитной проницаемости металла будем считать равным 1000. Для упрощения проверок удобно пронормировать параметры так, чтобы они изменялись от 0 до 1. В качестве точных будем использовать следующие значения параметров: 0.1,0.3,0.7,0.2,0.4,0.8.
Будем минимизировать отклонение значений магнитной индукции от значений полученных при точных параметрах в некотором наборе точек. Если обозначить через Bx,y,z,sx,sy,sz (р) - магнитную индукцию в точке р на наборе параметров Рисунок 4.16.: Поведение функционала в процессе минимизации на основе МКЭ в зависимости от метода оптимизации
Для нахождения минимума будем использовать метод градиентного спуска, метод сопряженных градиентов (Флетчера-Ривса) и квазиньютоновский метод BFGS (Бройдена-Флетчера-Гольдфарба-Шанно) [3; 4; 26; 62; 67].
Для аппроксимации решения в металле будем использовать метод конечных элементов. В воздухе будем использовать оба варианта аппроксимации. Поведение функционала при использовании в воздухе конечноэлементной аппроксимации изображено на рисунке 4.16. Из рисунка видно, что метод градиентного спуска сходиться значительно медленнее метода сопряженных градиентов, а метод BFGS до некоторого момента сходится примерно так же как и градиентный спуск, а затем скорость сходимости значительно возрастает, что соответствует теории. Отметим, что процесс минимизации с использованием метода конечных элементов является весьма чувствительным к сетке, а если точки, в которых сравнивается магнитная индукция, находятся в зависящей от параметров части сетки, оптими 0 5 10 — 15 20 25 30 35 40
Поведение функционала в процессе минимизации на основе МКЭ в зависимости от способа вычисления градиента зация, как правило, вообще не сходится из-за разрывности функционала (при ис Рисунок 4.19.: Сравнение поведения погрешности в зависимости от способа аппроксимации решения и метода оптимизации пользовании согласованных результант добиться устойчивой сходимости тоже не удается). На рисунке 4.17 изображено поведение функционала при использовании метода граничных элементов для аппроксимации решения в воздухе. Из рисунка видно, что поведение методов минимизации аналогично описанному выше.
Для вычисления градиента в описанных выше случаях использовалось численное дифференцирование. Если вся задача аппроксимируется с использованием метода конечных элементов, более выгодно использовать смешанный подход, изложенный в разделе 1.4. Сравнение сходимости процесса минимизации в зависимости от способа вычисления градиента приводится на рисунке 4.18. Из рисунка видно, что при использовании полуаналитического подхода сходимость несколько лучше, при этом время вычисления градиента с таким подходом может быть существенно меньше, особенно в случае нелинейной задачи.
Оптимизация геометрии ускорительных магнитов
Для построения базиса для пространства Я1 (П) требуется обеспечить непрерывность базисных функций. Базисные функции обычно получаются через суммы соответствующих локальных функций, определенных на элементах сетки. Процесс определения соответствия будем называть согласованием. Наиболее распространенным подходом к согласованию является ручной анализ базисных функций и основанная на этом анализе процедура согласования [18; 36; 38]. Такой подход оправдан, пока используется один базис, но при универсальной реализации метода конечных элементов предлагаемый подход с автоматическим анализом может позволить значительно упростить процедуры построения согласованных базисных функций.
Для построения базисных функций удобно ввести понятие мастер-элементов и определить на них локальные базисные функции. В качестве мастер-элементов обычно выбирают элементы максимально простой формы, например, единичный куб, что позволяет существенно упростить задание локальных базисных функций.
Базисные функции, с точки зрения согласования, можно разбить на 4 типа: 1. Вершинные - ассоциируются с вершиной сетки. Согласование обычно не представляет проблемы. 2. Реберные - ассоциируются с ребром сетки. Для согласования требуется учитывать направление ребра. 3. Граневые - ассоциируются с гранью сетки. Согласование требует учета ориентации грани и для функций высокого порядка является достаточно нетривиальной задачей.
4. Объемные, или внутренние, - ассоциируются с элементом сетки. Согласование для таких функций не требуется.
Для определения соответствия между мастер-элементом и элементом сетки необходимо сопоставить их локальную нумерацию вершин. В случае нерегулярной сетки соседние элементы могут быть сопоставлены с мастер-элементом достаточно произвольным образом, то есть локальные направления ребер на соседних элементах могут как совпадать, так и быть противоположны. Аналогичная ситуация и с ориентациями граней. Это накладывает определенные требования на построение базисных функций.
Минимальным требованием к базису для возможности его эффективного использования будем полагать неизменность пространства, порожденного базисом, при любой допустимой нумерации вершин. Это эквивалентно тому, что оператор перехода R может быть представлен в матричном виде: R( ф)=Щ (3.1) где г\) - вектор, составленный из базисных функций ф{, R - матрица оператора R. Будем рассматривать свойства локальных матриц мастер-элемента на примере матрицы жесткости. Теоретические выкладки для других матриц будут аналогичными. Матрица жесткости G имеет следующий вид:
Для произвольной матрицы R реализовать процедуру вычисления локальных матриц эффективным образом достаточно сложно. Наиболее удобным для последующей реализации является такой базис, где матрица R является матрицей перестановки. Это означает, что при перестановке вершин любая базисная функция переходит в какую-либо другую функцию этого же базиса. Таким свойством обладают, например, Лагранжевы базисы. Для построения согласованных базисных функций в этом случае достаточно при сопоставлении локальных и глобальных функций правильно учесть перестановку.
Однако построить иерархический базис, обладающий этим же свойством, невозможно, начиная с третьего порядка, поскольку функции нечетного порядка при смене направления ребра не могут отображаться в себя же. Для возможности эффективной реализации таких базисов необходимо учитывать, что при перестановке вершин базисные функции могут менять знак. Матрица оператора перехода R для таких базисов будет выглядеть следующим образом: R = PS, (3.4) где S - диагональная матрица, состоящая из 1 и —1, Р - матрица перестановки. Будем называть такие базисы перестановочными. Несмотря на то, что на практике используются и базисы не являющиеся перестановочными, универсальный и эффективный учет таких базисов достаточно сложен, поэтому в дальнейшем будем рассматривать только перестановочные базисы.
Рассмотрим вопрос автоматизированного анализа набора базисных функций на некотором мастер-элементе. Для этого, в первую очередь, необходимо определить, с чем ассоциирована каждая базисная функция. Для построения скалярных базисных функций из пространства Н1 () это можно сделать по следующему алгоритму: 1. Если базисная функция не равна нулю в какой-либо вершине элемента, то она ассоциируется с этой вершиной. 2. Если базисная функция равна нулю во всех вершинах элемента, но не равна нулю на каком-либо ребре, то она ассоциируется с этим ребром. 3. Если базисная функция равна нулю на всех ребрах элемента, но не равна нулю на какой-либо его грани, то она ассоциируется с этой гранью. 4. Если базисная функция равна нулю на всех гранях, то она ассоциируется со всем элементом.
Для других типов базисов можно построить аналогичную процедуру. Например, для трехмерных векторных функций пространства отличия будут заключаться в том, что вершинных функций среди них заведомо нет, а для проверки ре 59 берных и граневых нужно смотреть равенство нулю тангенциальной компоненты функции на соответствующих ребрах и гранях.
После определения ассоциаций требуется построить матрицы P и S. Эти матрицы будут иметь блочную структуру - они будут содержать отдельные блоки, отвечающие за вершинные функции, за реберные, за граневые, а также блок, отвечающий за объемные функции. Вершинные и объемный блоки будут единичными матрицами, поэтому интерес представляют реберные и граневые блоки. Для ребра возможно 2 направления, для треугольной грани - б ориентаций, а для четырехугольной - 8. Поэтому для эффективной реализации выгодно заранее, один раз для мастер-элемента, рассчитать все варианты блоков, а затем для элементов сетки выбирать требуемые блоки, исходя из нумерации вершин.
Определим для i-го ребра мастер-элемента преобразование системы координат R (е - переводящее это ребро в первое, и R-Де - переводящее это ребро в первое с изменением направления.
Таким образом можно вычислить матрицы перехода для инвертированного направления ребра. Если направление не инвертируется, то очевидно, что матрицы будут единичными. Для i-й грани мастер-элемента можно аналогично определить б преобразований, если грань треугольная, или 8 преобразований, если грань четырехугольная. Далее, если Rgce ( фі) ± R{fe (ijjj) = 0, (3.8) то фі переходит (возможно со сменой знака) при k-м варианте ориентации грани в ф3, т.е. Pj{ = 1, Sn = ТІ.
Таким образом можно вычислить блоки матриц перехода для всех вариантов ориентации ребер и граней. При реализации универсального конечноэлементного кода для различных базисов и мастер-элементов это позволяет значительно упростить разработку. Поскольку вычисления требуется проделать только один раз для базисных функций на мастер-элементе, дополнительные затраты времени на согласование не существенны.
Для вычисления локальных матриц в комплексе Quasar используется численное интегрирование. Для численного интегрирования треугольных и тетраэдральных элементов используются специальные квадратурные формулы, поскольку, несмотря на то, что квадратуры для них можно получить из одномерных, это является менее эффективным. Специальные квадратуры для треугольных элементов приводятся в статьях [37; 40; 84], для тетраэдральных элементов в [48; 86].