Диссертация (1149831), страница 5
Текст из файла (страница 5)
. . -ND, библиотека для C++, есть интерфейсы высокого уровня дляMatlab и для Python (!), поддерживает большое количество типов конечных элементов, включая экзотичные (например, X-FEM), практически любой размер-45ности, есть возможность программирования типовых задач с помощью готовых”кирпичиков”, избегая явной сборки линейной системы, отсутствуют встроенные средства создания сеток, можно пользоваться внешними, библиотека кроссплатформенна.LibMeshОписание:1-2-3D, библиотека для C++ с возможностью локальной адаптации сеток, параллельное решение линейных систем с помощью PETSc (MPI), поддерживаетбезматричные методы, выбор элементов шире, чем в deal.II.LifeVОписание:2D-3D, C++, основные области применения: гидродинамика, теплопроводность,перенос массы и взаимодействие жидкость–структура в пористых средах.OfeliОписание:2D-3D, библиотека C++, среди примеров есть задачи теплопроводости, решения уравнения Навье–Стокса для несжимаемой жидкости, теории упругости(2D и 3D), электромагнетизма, лицензия GPL.RheolefОписание:1-2-3D, библиотека для C++, код получается краток и приближен к математической записи задачи в слабой форме, автоматическая адаптация сеток для 2Dзадач.MODULEFОписание:Хорошо разработана библиотека для Fortran77.FEATFLOWОписание:462D-3D, библиотеки для Fortran, пакет ориентирован на решение уравненийНавье–Стокса для течения несжимаемых жидкостей, лицензия типа BSD.OpenFEMОписание:Хорошо разработана библиотека для Matlab и Scilab, но версия для свободногоScilab развивается менее активно.MelinaОписание:2D-3D, библиотека для Fortran, документирована исключительно на французском языке.FEAPpvОписание:библиотека для Fortran, для задач теории упругости и теплопроводности, распространяется бесплатно и в качестве приложения к книге, но лицензия неясна,является урезанной версией не свободной библиотеки FEAP.В диссертации в качестве ПО для решения задачи МКЭ выбранdeal.II [99, 100].
Пакет обладает свободной лицензией, наличием английской документации и примеров, позволяет решать 1-2-3D задачи.В deal.II не предусмотрен алгоритм создания сложной геометрии исследуемых систем и генерации сеток на ней, поэтому появилась необходимостьпоиска ПО для генерации сеток. Среди возможного ПО был выбран генераторсеток Gmsh, так как он единственный обладает свободной лицензией (GPL),графическим и командным интерфейсами, реализация на различных платформах(Windows, Linux, Mac OS X), так же данный пакет гибок и прост в использовании.Слабая формулировка задачиСлабая формулировка задачи может быть получена путем умножения уравнения (3.1) на тестовую функцию ∈ ′ (Ω), ∈ Ω ∈ 3 ,47 ′ = { () ∈ 2 (Ω)| () = 0, ∈ Γ } и интегрированием по граничной области Ω:∫︁−∫︁ · ∆ =Ω · .ΩИнтегрирование по частям дает:∫︁∫︁∇ · ∇ −· =ΩΓ∫︁ · .(3.2)ΩПреобразуем уравнение (3.2) к виду:∫︁∫︁∇ · ∇ =Ω∫︁ · +Ω · .ΓПодставив граничные условия Дирихле, получим:∫︁∫︁∇ · ∇ −Ω∫︁∇ · ∇ =Ω∫︁ · +Ω∫︁ · −Γ∇ · ∇ .ΩВоспользуемся решением ∈ 2 таким, что = − , с однородными граничными условиями = 0 ( ∈ Γ ).
Таким образом получим слабую формулировку уравнения Пуассона — найти такую функцию ∈ 2 , что выполняется∫︁∫︁∇ · ∇ =Ω∫︁ · +Γ∫︁ · +Ω∇ · ∇.(3.3)ΩДискретизация ГалеркинаДля дальнейшего решения задачи произведем дискретизацию наконечномерное пространство ⊂ 2 (Ω). Положим, что (1 , ..., ) являетсябазисом -мерного пространства и = ∩ 2 . Разложим функции (численное приближение ) и (численное приближение )по этому базису.48Таким образом, уравнение (3.3) приводит к следующим уравнениям:∑︁=1∫︁) · ∇ =(∇ +Ω∫︁ · −Γ∑︁∫︁∇ · ∇ . =1(3.4)ΩУравнение (3.4) представляет собой систему линейных алгебраических уравнений с постоянными коэффициентами: = ,где - симметричная, положительно определенная матрица.Для завершения дискретизации базис конечномерного пространства выбирается следующим образом: каждому узлу сопоставляется функция , которая равна 1 в и 0 для любого , ̸= .
Построение матрицы системы:двигаясь по узлам, строя базисные функции для каждого узла — определяютсялокальные матрицы, а затем, собирая локальные матрицы, получаем матрицусистемы = { }. Выражение для компонентов матрицы системы: =∑︁ ∫︁∈∇ ∇ где - множество узлов сетки.Построение сеткиДля построения сетки используется программный продукт Gmsh.Эта программа позволяет построить неструктурированную треугольную сеткуметодом Делоне, а затем разбивает её на четырехугольную. Для построенияфайла-сетки в Gmsh требуется выполнить следующие шаги:1) Задать геометрию системы:— задать точки основных геометрических узлов,49— соединить узлы кривыми (прямыми, секторами эллипсов, кривыми, определёнными пользователем).2) Выбрать алгоритм триангуляции, алгоритм разбиения элементов и др.
опции (в т.ч. формат файла-вывода).3) Выполнить команду для построения сетки и сохранения ее в файле.Алгоритм Делоне:1) Поместить узловые точки на углы области.2) Заключить всю область в границу.3) Триангулировать углы.4) Проверить, что триангуляция лежит в границе - 2).5) Поместить узловые точки в центры окружностей, описанных вокруг больших треугольников.6) Повторить из пункта 4), если Hmax(параметр, отвечающий за размер треугольника, наибольшая сторона) не достигнут.7) Удалить границы.Адаптивное улучшение сеткиДля адаптивного улучшения сетки используется реализация оценщика погрешности, разработанного Келли, Гаго, Зенкевичем и Бабушкой.Данныйметод аппроксимирует погрешность ячейки путем оценки интеграла скачка градиента через поверхность ячейки, функция оценки качества ячеек корректноприменима только для уравнения Пуассона с граничными условиями Дирихлеи ко-нормальной производной:1 с граничными условиями :1) Дирихле = ,2) Неймана · ⃗ · (∇ ) = .(︂)︂ 2+= − 250Функция оценки качества ячеекℓ2ℎ=24∫︁[ℎ]позволяет получить оценку для ячейки, показывая скачок градиента через сторону ячейки, где ℎ — наибольшая диагональ ячейки.
Если сторона ячейки принадлежит к типу границы Дирихле, то ее вклад не учитывается, если — Неймана, то оценивается следующим образом:∫︁[ −ℎ].После получения оценки для всех ячеек, улучшаются 30 процентов "самых плохих"из них.Метод сопряженных градиентов состоит в нахождении последовательности векторов ( ∈ ). Последовательность 1 ,2 ,..., является базисомточного решения:* = 1 · 1 + 2 · 2 + ...
+ · где = , = 1,2,..,. Точное решение доставляет минимум функционалу1() = − ,2который может быть переписан:∇() = − , ∇(* ) = 0.51Определим остаточный вектор линейной системы: = − = −∇( ), - направление градиента в . Новое направление спуска — +1 то же, чтои -сопряженный вектор с :+1 = − .
Алгоритм:∙ Инициализация:Пусть 0 будет начальным вектором,0 = − 0 ,1 = 0 ,0 =< 0 , 0 >.∙ Итерация: = 1, 2, . . . = · , = −1 / < , >, = −1 + , = −1 − , =< , >,+1 = +−1· .523.3Метод парных уравнений (МПУ)В данном параграфе для определения распределения электростатического потенциала в области, содержащей произвольное число соосных круговых диафрагм, используется метод парных уравнений (рис. 3.2).rR2 R3R1 R4Z1 Z2Z3 Z4zРис. 3.2: Схематическое изображение электронно-оптической системы счетырьмя круговыми диафрагмами.3.3.1Математическая модель системы соосных круговых диафрагмРаспределение потенциала должно удовлетворять уравнению Лапла-са с граничными условиями:(︂)︂⎧1 2⎪⎪+= 0,⎨ ⃒ 2⃒⎪⃒⎪(,±)= ± .⎩ ⃒ ≥ , - координаты диафрагм, с помощью которых моделируется фокусирующая косоугольная система.Так как (,) = − (, − ), то для решения задачи достаточнорассмотреть область ≥ 0.53Задача преобразуется к следующему виду:(︂)︂⎧1 2⎪⎪+= 0,⎪⎪ 2⎪⎨ (,0) = 0,⃒⎪⎪⃒⎪⎪⎪= .⎩ (, )⃒⃒ ≥ (3.5)Как решение граничной задачи (3.5), функция (,) может бытьпредставлена в виде разложения Ханкеля:+1 − (,) = +( − ) +−+1[︂∫︀∞ sh (+1 − )+ () +sh (+1 − )0]︂sh ( − )+1 () ×+sh (+1 − )×0 () ,(3.6) < < +1 .0 - функция Бесселя 0-го порядка.Итак, будем искать распределение потенциала в полуплоскости≥⃒⃒0, на границе которой потенциал удовлетворяет условию 0 = (,)⃒⃒= 0.=0Разобьем эту область на + 1 подобластей ≤ ≤ +1 .
На границах разделаподобластей = граничные условия на диафрагмах и условия непрерывности поля в отверстиях диафрагм имеют вид⎧⎪ 1 (,0) = 0; (, ) = +1 (, ), ≥ 0;⎪⎪⎪⎪⎨ (,∞) < ∞; (, ) = , ≥ ;⎪⎪⎪⎪ (,) +1 (,)⎪⎩=, = , < ;⃒⃒где +1 (,) = (,)⃒⃒, = 0, . ≤ ≤ +154Разложения Фурье-Ханкеля (3.6) удовлетворяют уравнению Лапласа, выполняется условие 0 () = 0 при = 0 и ограничивающее условие (,) < ∞ при → ∞, а также соблюдено требование непрерывности распределения потенциала.Граничные условия и непрерывность производной потенциала понормали к поверхности электродов (3.5) приводит к системе интегральныхуравнений:⎧ ∞∫︀⎪⎪ () 0 () = 0, ≤ < ∞;⎪⎪⎪0⎪[︃⎪⎪∞⎪∫︀1⎪⎪−1 () +−⎪⎪⎪sh(−)−10⎪⎪(︂⎪⎪ch ( − −1 )⎪⎪++⎨sh ( − −1)︂)ch (+1 − )⎪⎪+ () −⎪⎪⎪sh (+1 − )⎪]︃⎪⎪⎪1⎪⎪⎪−+1 () 0 () =⎪⎪sh(−)+1⎪⎪⎪⎪+1 − − −1⎪⎪=−, 0 ≤ < .⎩+1 − − −1Пусть⎧1⎪⎪()=−,⎪,−1⎪sh ( − −1 )⎪⎪⎪⎪ch ( − −1 )⎪⎪+⎨ , () =sh ( − −1 )ch (+1 − )⎪⎪+,⎪⎪sh (+1 − )⎪⎪⎪⎪+1 − − −1⎪⎪−.⎩ =+1 − − −1(3.7)55Тогда система (3.7) примет вид:⎧ ∫︀∞⎪⎪ () 0 () = 0, ≤ < ∞;⎪⎪⎪0⎪⎪∫︀∞ [︁⎪⎪⎨ −1, ()−1 () + , () () +0]︁⎪⎪++1, (+1 () 0 () = ,⎪⎪⎪⎪⎪⎪⎪⎩0≤<.(3.8)*Легко вычислить, что , () = 2(1 − ,()), где*,()1 (︁ exp(−(+1 − )+=−2sh (+1 − )exp(−( − −1 ) )︁+.sh ( − −1 )Очевидно, что выполняется условие:*lim ,() = 0.→∞Итак, система интегральных уравнений (3.8) может быть представлена как⎧ ∫︀∞⎪⎪ () 0 () = 0, ≤ < ∞;⎪⎪⎪0⎪⎪∫︀∞ [︁⎪⎪*⎨ −1, (−1 () + 2(1 − ,()) () +0]︁⎪⎪++1, ()+1 () 0 () = (),⎪⎪⎪⎪⎪⎪⎪⎩0≤<.(3.9)56Представим функции () в виде:∫︁ () = () sin .(3.10)0Рассмотрим разрывной интеграл Вебера∫︁∞0 () sin =0⎧⎪⎨ 0,0 ≤ < ;1⎪⎩ √, 0 < < .2 − 2Первые уравнения системы (3.9) с заменами (3.10) при использовании данногоразрывного интеграла тождественно выполняются:∫︀∞ ()0 () =0∫︀= () 0∫︀∞0 () sin ≡ 0, > .0Замены (3.10) для вторых уравнений системы (3.9) ведут к уравнениям:∫︀∞,−1 ()0 () 0−1 () sin +0+2+∫︀−1∫︀∞(1 −*,())0 () 0∞∫︀∫︀ () sin +0+1, ()0 () 0∫︀+1+1 () sin =0= (),0 ≤ < .(3.11)57В соответствии с методом решения парных интегральных уравнений [ []], уравнения (3.11) приводят к следующей системе уравнений:∫︀−1−1 () 0+2−2+∫︀0∫︀0∫︀+1 () () ∫︀∞0∫︀∞0∫︀∞,−1 ()1 () sin +1 () sin −*,()1 () sin +0+1 () 0∫︀∞,−1 ()1 () sin =0= (),где () =∫︀(3.12)0 ≤ < , .0Рассмотрим второй разрывной интеграл Вебера:∫︁∞1 () sin =0⎧⎪⎨ 0,0 ≤ < ;1⎪⎩ √, > .















