Содержание к диссертации
Введение
Глава 1. Вычислительные схемы моделирования двумерных электромагнитных полей. Построение двумерных конечноэлементных сеток. Двумерный препроцессор 22
1.1. Математические модели и вариационные постановки в двумерном случае 22
1.1.1. Модели двумерных нестационарных электромагнитных процессов... 23
1.1.2. Модели двумерных гармонических электромагнитных процессов 27
1.2. Математическое моделирование электромагнитного поля кабеля с корродирующей оболочкой 28
1.2.1. Постановка задачи 28
1.2.2. Вычисление синусоидальной и косинусоидальной компонент потенциала 30
1.2.3. Расчёт индуктивности электромагнитного поля кабеля с коррордирующей оболочкой при равномерном распределении тока по поверхностям жилы и оболочки 35
1.2.4. Расчёт индуктивности гармонического электромагнитного поля кабеля с корродирующей оболочкой 37
1.3. Пример решения задачи с использованием комбинированной сетки из треугольников и прямоугольников 42
1.3.1. Математическая модель процесса диссоциации двух зарядов 42
1.3.2. Конечноэлементная аппроксимация 43
1.4. Построение двумерных конечноэлементных сеток. Двумерный препроцессор 47
1.4.1. Описание расчётной области в двумерном препроцессоре 48
1.4.2. Проблемы выделения макроэлементов 60
1.4.3. Построение сетки на макроэлементе 65
1.5. Выводы з
Глава 2. Вычислительные схемы моделирования трехмерных нестационарных электромагнитных полей с использованием edge-элементов 78
2.1. Модель нестационарного электромагнитного поля в виде одного векторного уравнения 78
2.2. Вариационная и конечноэлементная постановка для модели с разрывным векторным потенциалом 81
2.3. Особенности конечноэлементной аппроксимации на edge-элементах... 86
2.3.1. Edge-элементы на параллелепипедах 86
2.3.2. Edge-элементы на тетраэдрах 87
2.3.3. Edge-элементы на призмах 91
2.4. Применение векторного МКЭ для анализа электромагнитного поля в согласованных плёночных СВЧ резисторах 93
2.4.1. Постановка задачи 94
2.4.2. Результаты численного моделирования 95
2.5. О решении гармонических задач векторным МКЭ 100
2.6. Выводы 102
Глава 3. Вычислительные схемы моделирования трехмерных нестационарных электромагнитных полей в технических устройствах с совместным использованием векторных и узловых элементов 104
3.1. Модель с совместным использованием вектор-потенциала электромагнитного поля и скалярного потенциала магнитного поля... 104
3.2. Вариационная постановка для модели с разрывным векторным и скалярным потенциалами 109
3.3. Конечноэлементная дискретизация для модели с разрывным векторным и скалярным потенциалами 115
3.4. Пример решения модельной задачи 118
3.5. Выводы 124
Глава 4. Вычислительные схемы моделирования трехмерных нестационарных электромагнитных полей с выделением главной части поля 125
4.1. Учёт нормального поля в системе уравнений Максвелла 125
4.2. Модель в виде уравнений для вектор-потенциала электромагнитного поля и скалярного потенциала электрического поля 127
4.3. Вариационная постановка и конечноэлементная аппроксимация 129
4.4. Учёт условий симметрии в постановке (4.27)-(4.30) 134
4.5. Модель и вариационная постановка с гармоническим источником поля 136
4.6. Выделение поля в модели с одним вектор-потенциалом 139
4.7. Выделение поля в модели с совместным использованием векторного и скалярного потенциалов магнитного поля 144
4.8. Использование вычислительных схем с выделением поля при решении задач геоэлектроразведки 147
4.9. Вычислительная схема численного моделирования электромагнитного зондирования Земли при наличии рельефа 151
4.10. Выводы 156
Глава 5. Вычислительные схемы моделирования нелинейных задач магнитостатики 158
5.1. Математическая модель на основе полного и неполного потенциала и конечноэлементная дискретизация , 159
5.2. Вычислительная схема с выделением главной части поля 164
5.3. Вычислительная схема для моделирования магнитного поля в конструкциях, имеющих хорошее двумерное приближение в декартовых координатах 169
5.4. Вычислительная схема для моделирования магитного поля в циклотронах 178
5.5. Выводы 183
Глава 6. Построение трёхмерных конечноэлементных сеток 185
6.1. Построение тетраэдральных сеток с использованием модификации метода тиражируемых сечений 185
6.1.1. Формирование конечноэлементной сетки по заданным сечениям 189
6.1.2. Построение тетраэдральной сетки по сетке из призм с треугольным основанием 190
6.1.3. Формирование информации о конечных элементах и об узлах сетки при проходе по сечениям 192
6.1.4. Краевые условия второго и третьего рода 193
6.2. Задание положения основных сечений в пространстве, описание их поверхностей и смещений узлов 194
6.3. Реализация рассмотренного подхода к построению трехмерных конечноэлементных сеток 197
6.4. Выводы 198
Глава 7. Алгоритмы, структуры данных и принципы построения конечноэлементного пакета TELMA 200
7.1. Особенности объектно-ориентированного конечноэлементного комплекса для моделирования электромагнитных полей 201
7.1.1. Необходимость выделения решаемой задачи как основной сущности201
7.1.2. Достоинства и недостатки объектно-ориентированного подхода к построению конечноэлементного комплекса 209
7.2. Объектно-ориентированная реализация алгоритмов работы с конечноэлементными СЛАУ в комплексе TELMA 211
7.2.1. Основные особенности конечноэлементных СЛАУ.. 212
7.2.2. Разреженный строчный формат хранения матрицы СЛАУ 213
7.2.3. Блочные форматы хранения. Разреженный строчно-блочный формат 216
7.3. Алгоритмы работы с конечноэлементными сетками 222
7.3.1. Построение портрета конечноэлементной матрицы .222
7.3.2. Нумерация глобальных базисных функций 226
7.3.3. Подсчёт и нумерация рёбер конечных элементов 229
7.3.4. Построение списков граней трёхмерных конечных элементов 231
7.4. Сборка глобальной матрицы из шаблонных локальных матриц 234
7.4.1. Явные и неявные схемы аппроксимации по времени 235
7.4.2. Шаблонные локальные матрицы для различных вариационных постановок 239
7.4.3. Сборка конечноэлементной СЛАУ для базового класса решения нестационарной электромагнитной задачи 242
7.5. Выводы 261
Глава 8. Вычислительные схемы для решения сопряженных задач 263
8.1. Постановка задачи и математическая модель .263
8.2. Вычислительные схемы 268
8.2.1. Вычислительная схема для совместного решения двух сопряженных нелинейных задач 268
8.2.2. Сравнение векторных элементов различных типов при решении задач с токами, текущими преимущественно в одном направлении 271
8.2.3. Вычислительная схема с несколькими двумерными моделями 276
8.3. Результаты численного моделирования. Сравнение с экспериментом.. 277
8.4. Выводы 281
Заключение 282
Список использованных источников
- Вычисление синусоидальной и косинусоидальной компонент потенциала
- Вариационная и конечноэлементная постановка для модели с разрывным векторным потенциалом
- Конечноэлементная дискретизация для модели с разрывным векторным и скалярным потенциалами
- Вычислительная схема численного моделирования электромагнитного зондирования Земли при наличии рельефа
Введение к работе
Актуальность темы. Для развития многих современных технологий очень важным является детальный анализ различных физических процессов, описываемых сложными математическими моделями. Изучение этих процессов невозможно без использования высокоточных методов численного моделирования. Для решения такого рода проблем к настоящему времени разработан большой объём программно-математического обеспечения, реализованного в широко известных многофункциональных программных комплексах, таких как ANSYS, NASTRAN, COSMOS, FLUX3D, OPERA3D, FEMLAB и т.д. Успешность этих комплексов определяется, в основном, тем, что в качестве базового метода моделирования в них используется метод конечных элементов (МКЭ). Тем не менее, на сегодняшний день и с помощью этих комплексов не удаётся получать результаты требуемого качества при решении многих очень важных практических задач.
К таким задачам можно отнести задачи моделирования трёхмерных электромагнитных полей в различных технических устройствах (установках индукционного нагрева и закаливания изделий, электролитических установках, электродвигателях и генераторах, ускорителях элементарных частиц и т.д.), а также задачи восстановления свойств объектов по данным электромагнитных зондирований (задачи дефектоскопии, задачи геоэлектроразведки и т.д.). Действительно, в указанных программных комплексах для моделирования трёхмерных полей заложены, в основном, достаточно простые модели (например, задачи электростатики, описываемые эллиптическим уравнением) и вычислительные схемы. Использование же таких программных комплексов для решения достаточно сложных трёхмерных задач, требующих разработки специализированных вычислительных схем, сторонним пользователем (т.е. без обращения к разработчикам) практически невозможно.
Например, в задачах геоэлектроразведки для восстановления проводимости среды по измеренным характеристикам электромагнитного поля необходимо получать решение соответствующих нестационарных трёхмерных задач моделирования электромагнитного поля с погрешностью до нескольких процентов в широком диапазоне времён во всех точках измерения, что не может быть обеспечено ни одним из перечисленных выше программных комплексов, использующих стандартные вычислительные схемы. Для решения таких задач необходимы специальные методы, позволяющие выделять поля влияния трёхмерных неоднородностей, и реализующие эти методы вычислительные схемы.
Большую сложность для моделирования представляют также трёхмерные задачи, в которых необходимо достаточно точно учитывать влияние вихревых токов на изучаемый физический процесс (например, индукционный нагрев или паразитные токи в электромеханических устройствах). Предлагаемые в настоящее время для решения таких задач подходы с использованием элементов векторного типа хотя и существенно расширяют класс решаемых задач, но имеют серьёзные проблемы при наличии в трёхмерной расчётной области непроводящих подобластей.
Необходимо отметить, что в России моделированию трёхмерных электромагнитных полей в сложных областях уделяется явно недостаточное внимание. Большинство современных работ в области численного моделирования электромагнитных полей направлено на решение либо одно- и двумерных задач (например, работы М.И.Эпова, В.С.Могилатова, Ю.А.Дашевского и др. ), либо достаточно простых трёхмерных задач (например, работы М.И.Эпова, Соловьёва С.А., Вабищевича П.Н., В.И.Дмитриева, М.Н. Бердичевского). При этом для решения некоторых (достаточно узких) классов трёхмерных задач разрабатываются быстрые вычислительные схемы, позволяющие приближать решение трёхмерной задачи решениями серии одно- или двумерных задач (например, работы М.И.Эпова, Ю.А.Дашевского). Программных же комплексов конечно-элементного моделирования электромагнитных полей, разработанных в России, вообще очень мало. Примерами могут служить программа ELCUT, разработанная в НПО ТОР (Санкт-Петербург) и предназначенный для решения двумерных нестационарных задач, а также программа MERMAID, разработанная А.Дубровиным в ИЯФ СО РАН и предназначенная для решения двумерных и трёхмерных задач магнитостатики на регулярных сетках. Эти пакеты являются достаточно узконаправленными, реализованные в этих пакетах вычислительные схемы не позволяют, например, моделировать трёхмерные нестационарные электромагнитные поля.
Таким образом, проблема построения эффективных вычислительных схем конечноэлементного моделирования, позволяющих получать численные решения задач моделирования электромагнитных полей в областях со сложной геометрией, и их гибкой настройки на особенности решаемой задачи до сих пор вызывает большой интерес у очень многих исследователей, занимающихся как вопросами разработки вычислительных методов, так и использующих эти методы при решении конкретных практических задач. Этим и определяется акту-. альность данной диссертационной работы.
Одной из важнейших при решении различных задач с использованием МКЭ является проблема построения конечноэлементных сеток. Значительное внимание, уделяемое исследователями этой проблеме, объясняется тем, что очень часто от того, насколько эффективно удается выполнить дискретизацию расчетной области, зависит сама возможность решения конкретной практической задачи с нужной точностью. Поэтому столь велик интерес исследователей к различным процедурам построения двумерных и трехмерных сеток. Этот интерес связан с тем, что именно использование неравномерных неструктурированных сеток часто позволяет при фиксированном числе узлов существенно уменьшить ошибку аппроксимации по сравнению с различными структурированными (например, квазипараллелепипеидальными) сетками. Повышение же точности численного расчета за счет уменьшения ошибок аппроксимации, связанных с дискретизацией исходной задачи, и минимизация вычислительных затрат за счет уменьшения размерности системы аппроксимирующих уравнений (а в МКЭ размерность этой системы напрямую зависит от числа узлов в конечноэлементной сетке) являются важнейшими проблемами при проведении любого численного исследования. Поэтому актуальной является разработка
вычислительных схем и алгоритмов, позволяющих использовать и строить как согласованные, так и несогласованные нерегулярные сетки, содержащие как однотипные, так и разнотипные элементы. Кроме того, при замене типа конечных элементов сама вычислительная схема практически не меняется, и поэтому желательно иметь программные инструменты, позволяющие достаточно гибко настраивать программный комплекс как на использование другой вычислительной схемы без изменения типа конечных элементов, так и на работу той же самой вычислительной схемы с другим типом элементов. Достаточно удобным для этих целей является применение объектно-ориентированного подхода при проектировании программного комплекса.
Развитие объектно-ориентированного программирования и появление достаточно надёжных и мощных компиляторов C++ закономерно привело к появлению объектно-ориентированных программных комплексов конечноэлемент-ного моделирования как в России, так и за её пределами. Однако, указывая очевидные преимущества объектно-ориентированных программных комплексов, такие как упрощение отладки и простота расширения пакета открытой структуры, большинство разработчиков решают, в основном, актуальную в задачах прочности проблему совместного использования разнотипных элементов, используя, фактически, одну вычислительную схему. При решении же трёхмерных задач электромагнетизма часто требуется одновременное использование нескольких вычислительных схем и вариационных постановок (например, при решении задач с выделением нормального поля одновременно решаются как минимум две связанные задачи: задача расчёта нормального поля и задача расчёта аномального поля, и для решения каждой из этих задач используется свой ' метод построения конечноэлементной аппроксимации).
В диссертационной работе много внимания уделено построению структур данных и алгоритмов, используемых при построении вычислительных схем ко-нечноэлементного моделирования электромагнитных полей, и анализу их эффективности при решении задач различных типов. Эффективность разработанных схем и их программной реализации будет продемонстрирована на примере решения сложных практических задач из различных областей науки и техники.
Таким образом, основной научной проблемой, пути решения которой рассматриваются в предлагаемой диссертационной работе, является проблема разработки и исследования эффективности методов конечноэлеменгного моделирования электромагнитных полей в сложных трехмерных областях.
Цель исследования состоит в разработке новых и повышении эффективности наиболее часто используемых перспективных методов численного моделирования трёхмерных электромагнитных процессов в сложных трёхмерных областях и программных средств, реализующих эти методы.
На защиту выносятся:
1) Методы конечноэлементного моделирования трехмерных нестационарных
электромагнитных полей с выделением главной части поля с использованием
узловых и векторных конечных элементов, а также при совместном их
использовании при исследовании электромагнитных процессов в
устройствах, состоящих из проводящих и непроводящих электрический ток конструктивных элементов.
Методы и средства описания двумерных и трёхмерных расчетных областей со сложными границами и алгоритмы построения двумерных и трёхмерных конечноэлементных сеток, реализованные в препроцессоре комплекса TELMA.
Методы конечноэлементного моделирования нелинейных магнитных полей в ускорителях заряженных частиц, использующие технологию поэтапного выделения главной части поля
Методы совместного решения нескольких взаимозависимых нелинейных задач электромагнетизма и теплообмена.
Объектно-ориентированная реализация библиотеки классов конечноэлементного моделирования в программном комплексе TELMA.
Научная новизна работы состоит в следующем.
Предложена конечноэлементная постановка с совместным использованием скалярных и векторных конечных элементов, основанная на вариационной постановке с использованием векторного потенциала в проводящих и скалярного в непроводящих подобластях. В отличие от постановок на основе только векторных конечных элементов, предложенная постановка при решении трёхмерных задач моделирования нестационарного электромагнитного поля без учёта токов смещения позволяет получать невырожденную систему линейных алгебраических уравнений даже при наличии в расчётной области непроводящих подобластей.
Разработаны новые средства описания сложной трехмерной геометрии и построения тетраэдральной конечноэлементной сетки с возможностями сгущения и разрежения её узлов. Эти средства позволяют при моделировании' электромагнитных полей в наиболее высокотехнологичных ускорителях заряженных частиц максимально точно учитывать конфигурацию фасок и шимм тем самым существенно улучшить качество получаемых результатов.
На основе конечноэлементных постановок с использованием векторных конечных элементов разработаны и реализованы вычислительные схемы с выделением главной части поля, позволяющие существенно повысить эффективность численного моделирования нестационарных электромагнитных полей в технических устройствах.
Впервые предложен и обоснован метод конечноэлементного моделирования нелинейных магнитных полей с выделением главной части поля. На его основе предложены и реализованы вычислительные схемы с поэтапным выделением главной части поля, позволяющие на порядок и более повысить точность численного моделирования нелинейных магнитных полей в ускорителях заряженных частиц.
Предложена и реализована объектно-ориентированная технология разработки программного комплекса конечноэлементного моделирования электромагнитных полей, обеспечивающая гибкость его настройки на используемые математические модели и вычислительные схемы с учётом особенностей решаемой задачи. На основе этой технологии были, например, решены
взаимозависимые нелинейные задачи электромагнетизма и теплообмена при моделировании электроконтактного нагрева криволинейного цилиндрического трубчатого изделия. Практическая ценность работы и реализация результатов. Разработанные методы и алгоритмы реализованы в программном комплексе TELMA и широко применялись для решения многих сложных практических задач из различных областей прикладных исследований: теплофизики, геофизики, электро- физики, а также при разработке и оптимизации различных электротехнических устройств. В диссертационной работе приводятся несколько примеров решения практических задач:
моделирование гармонического электромагнитного поля кабеля с корродирующей оболочкой;
расчёт характеристик процесса диссоциации двух зарядов с использованием комбинированной сетки;
анализ электромагнитного поля в согласованных плёночных СВЧ резисторах;
моделирование становления поля от вертикальной электрической линии в слоистых средах с трёхмерными объектами;
моделирование становления поля от петлевого источника в слоистых средах с трёхмерными объектами с учётом рельефа местности;
высокоточное моделирование магнитного поля в циклотроне
моделирование электроконтактного нагрева криволинейного цилиндрического трубчатого изделия.
Результаты диссертационной работы использовались при выполнении более чем 20 научно-исследовательских работ (как госбюджетных, так и хоздоговорных), из них в последние годы при выполнении тематических планов НИР НГТУ:
НГТУ. 1.4.99 «Исследование волновых моделей диагностики мест нарушения однородности в кабельных сетях среднего и высокого напряжения», 1999-2001;
НГТУ. 1.2.04 «Математическое моделирование электромагнитных процессов», 2004-2005;
а также хоздоговорных работ (за последние пять лет):
х/д ПМт-2-00 «Конечноэлементное моделирование геоэлектрических объектов» (2000г., ГФУП СНИИГГиМС);
х/д ПМт-3-00 «Расчёты электромагнитных полей при гальваническом возбуждении для сложных геоэлектрических условий Сибири» (2000г., ГФУП СНИИГГиМС);
х/д ПМт-1-01 «Решение прямой и обратной задачи геотермии в условиях седиментации» (2001г., ГП Дальинформгеоцентр);
' х/д ПМт-2-01 «Конечноэлементные расчёты трёхмерных магнитных полей» (2001г., НИУ ИЯФ СО РАН);
х/д ПМт-5-01 «Реализация учёта изменения палеоплотности и мощности слоев во времени при решении задачи геотермии в условиях седиментации» (2001г., ГПДальинформгеоцентр);
х/д ПМт-2-03 «Оценка глубинности и оптимизация основных параметров моноимпульсного электромагнитного зонда» (2003г., ГФУП СНИИГГиМС);
х/д ПМт-2-04 «Разработка и исследование технологии широкополосного импульсного зондирования» (2004г., ГФУП СНИИГГиМС);
х/д ПМт-2-05 «Проведение трехмерного конечноэлементного моделирования для технологии поиска золотосодержащих кварцевых жил методом вызванной поляризации» (2005г., ГФУП СНИИГГиМС);
х/д ПМт-3-05 «Расчёт параметров датчика ускорения» (2005г., ООО «Сенсор Текнолоджис»);
х/д ПМт-5-05 «Математическое моделирование электромагнитных полей при проведении аэроназемных поисково-оценочных исследований в условиях Восточной Сибири и Якутии» (2005г., ГФУП СНИИГГиМС);
' х/д ПМт-6-05 «Проведение региональных геофизических работ в зоне сочленения Сибирской платформы, Западно-Сибирской плиты и Енисей-Хатангского прогиба с целью подготовки новых зон нефтегазонакопления» (2005г., ГФУП СНИИГГиМС);
х/д ПМт-14-06 «конечноэлементные исследования трёхмерных магнитных
полей дипольных магнитов» (2006г., ИЯФ СО РАН);
Кроме того, исследования были частично поддержаны грантами РФФИ:
01-02-16932-а «Новый механизм поверхностной проводимости и появления зарядов в жидких диэлектриках»
03-02-16214-а «Экспериментальные исследования нового механизма появления зарядов на границе жидкого и твердого диэлектриков»
Достоверность результатов подтверждается как решением модельных задач и сравнением результатов решения некоторых частных задач с результатами других авторов, так и сравнением результатов численного моделирования с экспериментальными данными.
Личный вклад. Автором предложена и обоснована конечноэлементная постановка с совместным использованием узловых и векторных конечных элементов для модели электромагнитного поля с использованием разрывного векторного потенциала в проводящих и скалярного магнитного потенциала в непроводящих электрический ток конструктивных элементах, и на основе этой постановки разработана вычислительная схема с выделением главной части поля, описываемого двумерной (осесимметричной) задачей.
Автором предложены и реализованы в препроцессоре комплекса TELMA методы и средства описания двумерных и трёхмерных расчетных областей со сложными границами и алгоритмы построения двумерных и трёхмерных ко-нечноэлементных сеток
Автором разработаны и реализованы вычислительные схемы для решения задачи моделирования электромагнитных процессов в кабеле с корродирующей оболочкой и задачи моделирования процесса диссоциации для двух зарядов.
Автором предложены и разработаны методы конечноэлементного моделирования нелинейных магнитных полей в ускорителях заряженных, использующие технологию поэтапного выделения главной части поля
Автором разработаны и реализованы вычислительные схемы совместного решения нескольких взаимозависимых нелинейных задач электромагнетизма и теплообмена.
Автором выполнено численное моделирование трёхмерных нелинейных 'магнитных полей для циклотрона АІС-144, нестационарных электромагнитных полей с учётом рельефа местности при обработке данных зондирований на хребте «Безымянный», нестационарного трёхмерного электромагнитного поля согласованного плёночного СВЧ-резистора, электромагнитных и тепловых полей при решении задачи электроконтактного нагрева криволинейного цилиндрического трубчатого изделия.
Автором предложена объектно-ориентированная технология разработки программного комплекса конечноэлементного моделирования электромагнитных полей и реализована в виде библиотеки классов конечноэлементного моделирования в программном комплексе TELMA.
Автором выполнены все реализации рассмотренных в работе вычислительных схем и алгоритмов. Все исследования эффективности вычислительных схем и их реализаций, результаты которых приведены в диссертационной работе, также выполнены лично автором.
Апробация работы. Основные результаты работы были представлены и докладывались на: Всесоюзной конференции по методам численного решения многомерных нестационарных задач математической физики (Арзамас-16, 1991); Международной геофизической конференции и выставке по разведочной геофизике SEG-EAGO (Москва, 1993); Международном совещании-семинаре по механике реагирующих сред и экологии (Томск, 1994); First Asian Computational Fluid Dinamics Conference (Hong Kong, 1995); Particle Accelerator Conference and International Conference on High-Energy Accelerators (Dallas,1995); Международной конференции PaCT-95 (Санкт-Петербург, 1995); Международной геофизической конференции и выставке Санкт-Петербург'95; Международной конференции "неклассическая геоэлектрика" (Саратов, 1995); Международной конференции "Математические модели и численные методы механики сплошных сред" (Новосибирск, 1996); Научно-техническом совещании "Геофизические методы при разведке недр и экологических исследованиях" (Томск, 1996); Международной геофизической конференции "Электромагнитные исследования с контролируемыми источниками" (С.-Петербург, 1996); Третьем сибирском конгрессе по прикладной и индустриальной математике, посвященном памяти С.Л.Соболева, 1998г; Международной геофизической конференции «300 лет горно-геологической службе России», Санкт-Петербург, ,2000г.; 63 и 65 Международной конференции EAGE Conference & Technical Exhibition, (Amsterdam, 11-15 June 2001, Stavanger, Norway, 2-5 June 2003); международной конференции по вычислительной математике МКВМ-2004; Шестой международной конференции «On UNCONVENTIONAL ELECTROMECHANICAL AND ELECTRIACAL SYSTEMS», 2004г.; Международной конференции
«International Symposium On Heating by Electromagnetic Sources», Padua, June 22-25, 2004г.; Третьей, четвёртой, седьмой и восьмой международной научно-технической конференции "Актуальные проблемы электронного приборостроения" АПЭП-96, АПЭП-98, АПЭП-2004, АПЭП-2006; Третьей и восьмой Российско-Корейской международной конференции (KORUS'99, KORUS-2004); а также на семинарах ВЦ СО РАН, МВТ СО РАН, ИТПМ СО РАН, ИЯФ СО РАН, ИВМиМГ СО РАН (г.Новосибирск), ОИЯИ (г.Дубна). Результаты автора включались в отчеты по НИР НГТУ, заключительные отчеты СНИИГ-ГиМС.
Публикации. Результаты исследований изложены в 56 печатных работах (32 статьях и 24 трудах российских и международных конференций), среди которых 20 статей опубликованы в изданиях, включенных в перечень ВАК ведущих рецензируемых научных журналов и изданий, выпускаемых в Российской Федерации, в которых должны быть опубликованы основные научные результаты диссертаций на соискание ученой степени доктора наук.
Структура работы. Диссертация состоит из введения, восьми глав, заключения, списка использованных источников (269 наименований) и приложения. Работа изложена на 319 страницах, включая 47 рисунков и 5 таблиц.
Вычисление синусоидальной и косинусоидальной компонент потенциала
Однако на практике это обычно приводит к заметным (1-2%) ошибкам в результатах моделирования (на достаточно поздних временах). Это связано с тем, что, задавая в узлах значение решения, мы фактически задаём веса конечноэлементной аппроксимации. Но конечноэлементное решение, получаемое при решении краевой задачи, является ближайшим к точному на сетке в интегральном смысле, а не в узлах, и поэтому эти веса, естественно, не совпадают с весами конечноэлементного решения. Вычислительные эксперименты показывают, что даже на достаточно грубой сетке, если в качестве начального приближения взять не значения в узлах (1.14), а решить методом конечных элементов на той же сетке краевую задачу для уравнения — div -grad + T = JP (116) \x,r с б-функцией в правой части, то погрешность моделирования на поздних временах (на которых, собственно, и требуется получать решение) окажется меньше. 7.7.2. Модели двумерных гармонических электромагнитных процессов В случае, когда в уравнении вида (1.9) правая часть зависит от времени гармонически, т.е. J{t) = Jssin(w) + Jccos(io), (1.17) решение уравнения (1.9) может быть найдено в виде u{t) — u5sin(w) + ис cos (cot), (1-18) где us и ис определяются из решения системы уравнений -div(\gradu5) + ca/ + стишс -euV = Js, (1.19) - div (X grad uc) + cxu - стиж5 - єш V = Jc. (1.20) Обратим внимание на то, что оба уравнения (1.19) и (1.20) фактически при конечноэлементной аппроксимации дадут те же самые интегралы по конечным элементам, что и уравнение (1.11) (см. [133]), чем мы воспользуемся в дальнейшем при построении универсального алгоритма сборки в программном комплексе.
Заметим, что для большинства реальных задач моделирования электромагнитных процессов универсальная технология решения краевых задач для уравнений вида (1.11) или (1.19) и (1.20) не является самоцелью, а служит только частью общей вычислительной схемы решения поставленной задачи, что требует достаточной гибкости разрабатываемых универсальных алгоритмов. Примером такой задачи является задача моделирования электромагнитного поля с целью разработки методики определения волнового сопротивления участка кабельной линии при коррозии металлической оболочки [101], которая была решена с использованием разработанных автором вычислительных схем.
Скрытые дефекты в силовых кабельных линиях среднего и высокого напряжения в ряде случаев могут привести к внезапному выходу их из строя и, соответственно, к существенному ущербу потребителей электроэнергии. К таким скрытым дефектам относится корро место коррозии оболочки кабеля
Эскиз участка кабельной линии с корродирующей металлической оболочкой приведен на рис. 1.1. При допущении об однородности зоны корродирования по всей длине поврежденный участок кабельной линии может моделироваться в виде однородной длинной линии. При этом допущении наиболее информативным параметром, позволяющим идентифицировать коррозию металлической оболочки и в ряде случаев даже определить её размеры, является погонная индуктивность поврежденного участка линии. Эта индуктивность определяется с помощью численного расчета электромагнитного поля. Поскольку проводимость грунта, в котором проложен кабель, существенно меньше проводимости металлической оболочки, можно считать, что кабель окружает диэлектрик.
Будем считать, что зона корродирования однородна по всей длине поврежденного участка. Тогда при достаточно большой длине этого участка можно считать, что электромагнитное поле описывается системой уравнений которые являются частным случаем системы (1.19)-(1.20) в декартовой системе координат при пренебрежимо малых токах смещения. При этом на удалённой границе расчетной области заданы краевые условия: -0. (1.23) Гі 0, А Расчётная модель приведена на рис. 1.2. Будем решать задачу при следующих исходных данных: Щ = 0.01м, Rx = 0.016 м, jRg = 0.01625 м. Будем считать, что удельная проводимость зо жилы и оболочки совпадает и равна 5.6-10 См/м, магнитные проницаемости жилы и оболочки совпадают с магнитной проницаемостью вакуума. Будем также считать, что кабель находится в непроводящей среде, а частота источника поля равна 10 Гц.
Фрагмент конечноэлементной сетки для расчета поля показан на рис.1.3. Для расчёта электромагнитного поля в кабеле необходимо задать сторонние токи в жиле и оболочке (токи, обусловленные напряжением). Определить их величину можно следующим образом. где 5Ж - площадь жилы, и Jz =0. Решив полученную краевую задачу, проинтегрируем решение по подобластям расчетной области, чтобы рассчитать наведенные токи. Обозначим STS ж ж =i af где 1AG, 4f -решение краевой задачи (1.21)-(1.23) с источником (1.24), Ож - подобласть расчётной области, соответствующая жиле, Поб - подобласть расчётной области, соответствующая оболочке. Верхний левый индекс S (или С) у токов означает, что токи получены при синусоидальном (или косину-соидальном) стороннем токе, а левый нижний индекс ж (или об) указывает, что сторонний ток задан в жиле (или в оболочке). Правый верхний индекс S (или С) указывает, что ток синусоидальный (или косинусоидальный), правый нижний индекс {ж или об) обозначает область интегрирования (жила или оболочка).
Вариационная и конечноэлементная постановка для модели с разрывным векторным потенциалом
При равномерном распределении постоянного тока по поверхностям жилы и оболочки существует аналитическое решение, позволяющее рассчитать величину относительной индуктивности в зависимости от величины угла поврежденного участка G [101]. Рассчитаем эту величину численно и сопоставим с результатами аналитического расчёта.
Электромагнитное поле в области, изображенной на рис. 1.2, при равномерном распределении постоянного тока по поверхностям жилы и оболочки рассчитывалось на треугольной сетке с 25100 узлов и 49924 треугольников.
Относительная индуктивность кабеля при величине угла поврежденного участка 6 вычислялась по формуле / BldVt Bfan Me) = J-7Tr О-38) п где BQ - модуль вектора магнитной индукции, рассчитанный при величине угла поврежденного участка G, В0 - модуль вектора магнитной индукции, рассчитанный при отсутствии повреждения.
На рис. 1.4 приведен график зависимости относительной индуктивности от угла поврежденного участка. Различия в значениях индуктивности, вычисленных по численному и аналитическому решениям, нигде не превосходили 0.1%.
Зависимость относительной индуктивности от угла поврежденного участка при равномерном распределении тока по поверхности жилы и оболочки 1.2.4. Расчёт индуктивности гармонического электромагнитного поля кабеля с корродирующей оболочкой
Рассчитаем теперь относительную индуктивность кабеля при условии, что реальные токи, наводимые в жиле и оболочке, могут быть распределены неравномерно по поверхностям жилы и оболочки. Находить эти токи будем с помощью описанной выше методики моделирования.
Относительной индуктивностью подобласти П при величине угла поврежденного участка 0 будем называть величину где BQ И BQ - соответственно синус-компонента и косинус-компонента модуля вектора магнитной индукции, рассчитанный при величине угла повреж-денного участка G, В0 иВ0 - соответственно синус-компонента и косинус-компонента модуля вектора магнитной индукции, рассчитанный при отсут-ствии повреждения (т.е. при 0 = 0). Пусть Д , AQ - потенциал, рассчитанный по формулам (1.37) для задачи с величиной угла поврежденного участка G. Тогда
На рис. 1.5 приведен график зависимости относительной индуктивности всей расчетной области от угла поврежденного участка. Отметим, что эта зависимость существенно отличается от зависимости, соответствующей постоянному распределению плотности тока по поверхности жилы и оболочки, приведённой на рис. 1.4. Покажем, что распределение токов по поверхности жилы и оболочки действительно неравномерно по окружности. Для этого разобьем подобласти расчётной области, соответствующие оболочке и поверхности жилы, на участки по 9 градусов, и проинтегрируем по этим участкам модуль суммарных токов. Полученные величины пронормируем на значения таких же интегралов по тем же подобластям для токов при отсутствии повреждения. На рис. 1.6 и рис. 1.7 приведены зависимости пронормированных модулей полного тока от величины угла, отсчитываемого от вертикальной оси по часовой стрелке (см. рис. 1.2). Заметим, что увеличение поврежденного участка оболочки кабеля приводит к нарастанию неравномерности распределения тока не только по поверхности неповрежденного участка оболочки, но и по поверхности жилы.
Так, при повреждении оболочки кабеля наполовину (коррозия 180) плотность тока на поверхности жилы, находящейся напротив центра неповрежденного участка оболочки (т.е. при ц? = 180), более чем в три раза превосходит плотность тока на поверхности жилы, находящейся напротив центра поврежденного участка оболочки (т.е. при ip = 0). Относительная же индуктивность для этого случая (т.е. при повреждении оболочки кабеля наполовину) Ln (90) = 1.9 при равномерном распределении тока по поверхности жилы и оболочки и Ln(90) = 1.5 при неравномерном (вычисленном с учётом взаимоиндукции) распределении тока в жиле и оболочке.
Необходимо отметить, что относительная индуктивность кабеля практически неизменна в достаточно широком диапазоне частот. Так, для частоты 1МГц относительная индуктивность при любом угле поврежденного участка отличается от относительной индуктивности для частоты 100МГц не более чем на 0.2%, а для частоты 10МГц - не более чем на 0.05%.
Конечноэлементная дискретизация для модели с разрывным векторным и скалярным потенциалами
Будем перемещать фронт генерации сетки по построенным линиям, заполняя каждую получившуюся полосу треугольниками. При заполнении полосы треугольниками каждый следующий треугольник получается присоединением новой вершины с одной из сторон полосы к крайнему ребру уже построенных треугольников (при этом сторона, с которой берётся вершина, выбирается так, чтобы получившийся треугольники имел менее острые углы - это делается сравнением длин соответствующих отрезков, как это будет показано на примере).
Процесс построения сетки проиллюстрирован на рис. 1.11. Занумеруем узлы сетки так, как это показано на рис. 1.106. Поясним процесс заполнения треугольниками первой полосы. Процесс начинается с ребра 1-2. Содержа 69 щий это ребро треугольник может быть построен двумя способами: либо подсоединением этого ребра к вершине 5 с левой стороны полосы (образуется треугольник 1-2-5), либо к вершине 6 с правой стороны (образуется треугольник 1-2-6). Из этих двух треугольников менее острыми углами обладает тот, который образуется более коротким ребром из рёбер 1-6 и 2-5, т.е. 1-2-6 (рис. 1.11а). Процесс продолжается для ребра 1-6. Из двух возможных способов подсоединения - ребром 6-5 или ребром 1-12, выбирается подсоединение более коротким ребром 6-5, как это показано на рис. 1.116. Теперь фронтальным ребром становится ребро 6-5 и т.д. Результат построения сетки показан на рис. 1.1 Із.
Пусть на стороне у=0 задано уо внутренних узлов, на стороне у=\ задано vj/i узлов, на стороне х=0 задано %о внутренних узлов и на стороне х=1 задано ,i внутренних узлов. Пусть, для определенности, 6,0 6,l5 ф0 і\)г и выполняется либо условие: — —, либо, при — = —, фр Основная ,і Фі і Фі идея алгоритма заключается в следующем. От единичного квадрата последовательно отсекаются трапеции с основаниями, являющимися частями сторон у=0 и у=1 рассматриваемого квадрата, так, чтобы боковые стороны этих трапеций содержали одинаковое число узлов. Для генерации сетки в трапециях можно использовать алгоритм, описанный выше. Рассмотрим отсечение трапеций более подробно. 6 5
Попытаемся построить сетку с почти линейно меняющимся вдоль оси X распределением узлов. Для этого будем отсекать от единичного квадрата трапеции так, чтобы их боковые стороны были наиболее близки к тем параллельным оси Y линиям, на которых количество узлов должно измениться. При этом вершины трапеций должны совпадать с узлами на границах обрабатываемого квадрата.
Пусть JCJ - координата линии, при пересечении которой количество узлов должно измениться с / на (г+1). Тогда номер узла на стороне у=0, который должен использоваться для построения трапеции, можно определить по формуле А4=Мф0+1) + 0.5 Аналогично определяется номер узла Ц на стороне у=\: fcj = \Х{ (ф1 + 1) + 0.5 (1.58) При этом значение Х{ можно вычислить как тк±! ti-e.o+1" Однако, если для построения вершин отсекаемой трапеции использовать только эти формулы, можно не добиться желаемого результата, поскольку отсутствует запрет на получение треугольников. Кроме того, иногда возможно получение сетки с резким изменением шага. Например, если на одном основании трапеции окажется один внутренний узел, а на другом - ни одного, то шаг изменится вдвое. Поэтому необходимо произвести корректировку вычисленных значений к$ и fc[, учитывая значения fcg и к для предыдущего отсечения. Эта корректировка х\ может быть выполнена по формуле: Выполним следующие корректировки: 1. Если &о &о , положим &Q = &о , изменим Х[, используя (1.59), и пересчитаем Ц по формуле (1.58). 2. Если Q А 1"1, положим = к[ 1. 3. Если кг0 = klQ1 + 1 и kl = kl l + 2, положим к$ = к$ +2, изменим х\, используя (1.59), и пересчитаем fc{ по формуле (1.58). 4. Если fcj = fcjl +1 и fcj к{"1 + 3, положим fcj = к%l + 2, 5. Если после всех корректировок &Q = Ц" или & = fc" , то текущая трапеция не строится (поскольку она вырожденная). Сделанные корректировки могли сместить боковую сторону трапеции в область, где количество внутренних узлов на ней должно уже увеличиться. Поэтому необходимо изменить значение / по формуле і = max fa (І, - І0 + 1) + І0 - 1, і). (1.60)
Для того, чтобы избавиться от резких переходов шага, можно использовать для генерации узлов внутри трапеции не полученное по формуле (1.60) число узлов і на боковой стороне, а число т\ определяемое как: Пусть задан многоугольник, для четырех вершин которого А = (Ах,Ау), В = (Вх,Ву), С — (Сх,Су), D == (Dx,Dy) задано соответствие вершинам четырехугольника, а для вершин четырехугольника - соответствие вершинам единичного квадрата. Будем считать, что на каждой стороне этого единичного квадрата задано число внутренних узлов, равное количеству вершин исходного многоугольника между двумя его вершинами, соответствующими вершинам рассматриваемой стороны единичного квадрата. Для генерации узлов внутри многоугольника можно построить отображение F, переводящее узлы, сгенерированные в единичном квадрате, на внутренность многоугольника. Будем строить отображение F так, чтобы прообразом вершины А многоугольника была точка (0,0), то есть F l (А) = (0,0); аналогично для остальных вершин многоугольника: F l (В\ = (0,1), F l (С) = (1,1), F-1 (3) = (1,0).
Вычислительная схема численного моделирования электромагнитного зондирования Земли при наличии рельефа
В качестве примера использования векторного МКЭ рассмотрим задачу анализа электромагнитного поля в согласованных плёночных СВЧ резисторах [150].
В СВЧ диапазоне для построения мощных широкополосных нагрузок применяют плёночные резисторы сосредоточенного или распределённого типа. В плёночных резисторах мощностью 100Вт и более на невысоких частотах (порядка сотен мегагерц) ток, протекающий по резистивной плёнке, практически не имеет пространственной дифференциации, из-за относительно невысокой удельной проводимости плёнки и пренебрежимо малых токов смещения. Поэтому при определении входного импеданса и источников тепловыделения для моделирования температурных полей достаточно использо 94 вать относительно простые двумерные математические модели электромагнитного поля. Для частот же порядка 1 ГГц и выше ток в плёнке изменяется не только (и даже не столько) в поперечном сечении, но и вдоль неё (из-за существенного влияния токов смещения). В этом случае достаточно точные оценки характеристик резистора можно получить только на основе трёхмерного моделирования электромагнитного поля.
Математическое моделирование электромагнитных процессов будем проводить для конструкции, состоящей из резистивной плёнки размером 6 х 8 х 0.006 мм3 с удельной проводимостью о = 4.17 103См/м и двух медных контактов размерами 6 х 10 х 0.006 мм3 и 6 х 1 х 0.006 мм3, нанесённых на диэлектрическую подложку размером 10 х 19 х 4 мм3 с диэлектрической проницаемостью є = 6є0 (є0 - диэлектрическая проницаемость вакуума). Сторонний ток, имитирующий генератор, задан в отдельной подобласти размером 1 х 0.1 х 4 мм , соединяющей медный контакт с идеально проводящим основанием. Второй медный контакт выходит на идеальный проводник, соединяющий резистор с идеально проводящим основанием. Схематический рисунок конструкции приведен на рис.2.3. Вся конструкция закрыта идеально проводящим экраном, удалённым от резистора на 60 -г- 70мм. Поскольку задача моделирования электромагнитного поля имеет плоскость симметрии х = 0, разрезающую конструкцию на две одинаковые части, в расчётную областью может быть включена только половина описанной конструкции. На этой плоскости в качестве условия симметрии необходимо задать однородное краевое условие (2.9).
На тех частях SE границы S, которые соответствуют поверхности идеального проводника, должно быть задано краевое условие Ёхп = 0, (2.39) — т.е. касательные составляющие Е на поверхности идеального проводника должны равняться нулю. Соответственно, на SE должно выполняться анало —» — гичное уравнение и для вектор-потенциала А и пробных вектор-функций Y, что соответствует однородному краевому условию (2.8). Численные расчёты проводились по изложенной выше вычислительной схеме для гармонических сторонних токов с частотами в диапазоне от 200МГц до 5ГГц.
Результаты, полученные для частот до 500МГц, практически не отличались от результатов двумерного моделирования, когда считалось, что токи текут только параллельно продольной оси резистора, а магнитное поле имеет только перпендикулярные к этой оси составляющие. На более же высоких частотах, когда с проводников стекают существенные токи смещения, картина электромагнитного поля может кардинально отличаться от полученной с помощью двумерного моделирования. Мы приведём результаты численного моделирования электромагнитного поля для частоты 5ГГц. Величина стороннего тока была подобрана (для рассматриваемой линейной задачи это просто множитель) так, чтобы выделяемая на резистивной плёнке тепловая мощность составляла 200Вт.
На рис.2.4 представлено распределение вектора плотности тока в продольном сечении резистивной плёнки, перпендикулярном плоскости симметрии, в моменты времени, когда сторонний ток максимален (а) и равен нулю (б). Из этого рисунка видно, что поле токов очень дифференцированно не только поперек, но и (и даже более существенно) вдоль плёнки. Это объясняется существенным влиянием токов смещения, распределение вектора плотности которых в плоскости симметрии х = 0 приведено на рис. 2.5 (рис. 2.5а соответствует моменту времени, когда сторонний ток максимален, а рис. 2.56 - когда сторонний ток равен нулю). Распределение вектора поверхностной плотности тока Jn0B в те же моменты времени на идеально проводящем основании представлено на рис.2.6. Поверхностная плотность тока на идеальном проводнике определяется соотношением:
