4.2 Препроцессор

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

В соответствии с особенностями МКЭ на этапе препроцессора пользователь должен выполнить ряд основных шагов:

– задание уравнения в частных производных;

– ввод геометрии области, разбитой на ряд простых подобластей с учетом особенностей решаемой задачи;

– задание граничных и других условий, неоднородностей среды;

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

– сохранение данных на диске.

Препроцессор состоит из нескольких функциональных блоков, которые приведены на рис. 4.1. Примерная последовательность действий при работе с препроцессором показана на рис. 4.2.

Блок построения геометрии области. Постановка краевой задачи всегда подразумевает задание пространственной области, на которой решается дифференциальное уравнение

Рис. 4.1. Обмен данными в препроцессоре.

Рис. 4.2. Последовательность действий при работе с препроцессором.

вместе с граничными условиями. Затем в соответствии с МКЭ эта область должна разбиваться на конечные элементы без пропусков и перекрытий. Поэтому построение расчетной области задачи должно учитывать используемые алгоритмы такого разбиения (триангуляции). В рассматриваемом пакете реализован способ триангуляции с предварительным разбиением области на сравнительно крупные и сравнительно простые подобласти и последующим их разбиением на конечные элементы. Такой подход облегчает контроль качества конечно-элементной сетки, управление плотностью распределения узлов, а также сильно упрощает алгоритм дискретизации.

Принципы построения области, заложенные в препроцессор, следующие. Имеется три базовых геометрических объекта: точка (вспомогательный узел), линия и зона. Вспомогательные узлы задаются в виде координат (пар чисел). Каждому введенному узлу ставится в соответствие номер, который отображается на экране вместе с соответствующей точкой. По заданным точкам определяется линия. Это может быть либо прямая, либо кривая 2-го порядка, либо дуга окружности. Чтобы провести прямую, точнее, отрезок, соединяющий две точки, пользователь должен ввести два номера уже существующих узлов или координаты начала и конца отрезка. Через три заданные точки (x₁,y₁), (x₂,y₂), (x₃,y₃) проводится дуга кривой, параметрическое уравнение которой записывается в виде

(4.1)

где параметр t лежит в интервале [–1, 1]. При этом пользователем вводятся номера узлов, совпадающих с этими точками, либо их координаты. Ввод дуги окружности осуществляется заданием начального и конечного узла, а также центра окружности. Введенная линия отображается на экране вместе со своим уникальным номером.

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

В программе можно задавать зоны двух типов – четырехугольные и зоны Делоне. Такое деление обусловлено двумя различными методами триангуляции, имеющимися в препроцессоре. Для разбиения зон первого типа используется блочный метод и метод Делоне, для зон второго типа – только метод Делоне. Предпочтение отдается четырехугольным зонам: они обеспечивают более равномерное и предсказуемое разбиение.

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

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

Пусть, например, расчетная область задачи представляет со­бой прямоугольник с круглым отверстием в середине (рис. 4.3).

Рис. 4.3. Расчетная область задачи

Из соображений симметрии достаточно построить четверть области и разбить ее на зоны (рис. 4.4).

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

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

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

Рис. 4.4. Ввод геометрии области.

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

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

Блок задания параметров задачи. Данный блок позволяет определить тип дифференциального уравнения, которому удовлетворяет искомая функция; стационарная или нестационарная задача; двумерная или осесимметричная.

В общем случае дифференциальное уравнение имеет вид:

,

где – неизвестная функция, вид функций G, A и F определяется пользователем. Данное уравнение имеет достаточно общий вид. В частном случае, когда G=F=0, =const получим уравнение Лапласа:

;

если G=0, , =const, то уравнение Пуассона:

;

если G=0, =const, , то уравнение Гельмгольца:

;

если G=0, =const, , то уравнение Лондонов:

;

если , , то уравнение теплопроводности:

;

если , , то уравнение диффузии:

.

Кроме этого, в данном блоке осуществляется ввод граничных и других условий. На линии Г могут быть заданы либо граничное условие Дирихле (условие 1-го рода):

,

либо условие Неймана (2-го рода):

,

либо смешанное (условие 3-го рода):

.

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

u₊ – u_– = I,

где u₊ и u_– – значения потенциала по разные стороны от линии разреза. Величина скачка потенциала численно равна току, протекающему по сверхпроводнику. Выбирать линию разреза следует таким образом, чтобы любой замкнутый контур, охватывающий сверхпроводник с током, обязательно пересекал разрез. Кроме величины скачка, требуется еще указать направление обхода, при котором скачок потенциала положителен (различные направления обхода соответствуют разным направлениям тока, текущего по сверхпроводнику).

Кроме того, доступны следующие дополнительные типы условий на линиях:

• условие постоянства неизвестной функции:

u=const;

• интегральное условие:

;

• условие заданного потока через линию разреза Г₁ для многосвязной области:

(знак «+» или «–» определяется направлением нормали к поверхности Г₁).

В ряде случаев математическая формулировка допускает существование решения для потенциала с точность до произвольной постоянной (например, в задачах с условием «разрез»). Для численного метода, каким является МКЭ, такая неоднозначность неприемлема. Поэтому для выделения решения необходимо задать фиксированное значение в каком-либо узле конечно-элементной сетки:

u_k=U.

Для нестационарных задач необходимо также задать начальное условие – значения функции в момент времени t = 0:

.

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

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

– константа f = C;

– полином 1 порядка ( f = Ax + By + C или f = Au + B) или более высокого порядка;

– экспонента f = Aexp(Bx + Cy + D) или f = Aexp(Bu + C);

– тригонометрические функции (синус или косинус) f = Asin(Bx + Cy + D) или f = Acos(Bu + C);

– степенная (f = Ax^B +Cy^D или f = Au^B).

Эти элементарные функции являются функциями только одной переменной (либо только , либо только u, и т. д.). Они комбинироваться между собой операциями сложения (f₁ + f₂), умножения (f₁ * f₂) или деления (f₁ / f₂).

Блок визуализации. Функции блока визуализации обеспечивают отображение вводимой области задачи в удобном для пользователя виде. Предусмотрены:

– сдвиг изображения и изменение его масштаба;

– разрешение/запрещение отображения номеров объектов;

– отображение краевых и других условий,

– отображение конечно-элементной сетки, номеров ее узлов и элементов, и др.

– отображение прямоугольной координатной сетки.

Блок триангуляции подобластей. Генерация конечно-элементной сетки является наиболее важным этапом работы препроцессора, так как слишком редкая или некачественная сетка может привести к большим погрешностям решения задачи, а слишком частая – к неоправданным затратам машинного времени на этапе вычислений. Разбиение должно быть проведено без пропусков и перекрытий, при этом форма конечных элементов должна не очень сильно должна отличаться от правильной. При проведении разбиения следует учитывать априорную информацию о поведении неизвестной функции u: там, где предполагается быстрое изменение u (например, вблизи неоднородностей, острых углов и т. п.), сетка должна быть более густой, там, где изменение u должно быть медленным – более редкой. На самом деле эта информация должна учитываться уже на этапе предварительного разбиения области на зоны.

В программе FEMPDESolver для решения двумерных и осесимметричных задач могут использоваться конечные элементы следующих типов:

– лагранжев треугольный элемент 1-го порядка (3 узла);

– изопараметрический треугольный элемент 2-го порядка (6 узлов);

– серендипов четырехугольный элемент 1-го порядка (4 узла);

– изопараметрический четырехугольный элемент 2-го порядка (8 узлов);

– треугольный эрмитов элемент (4 узла, 10 степеней свободы);

– бесконечный элемент 1-го порядка (4 узла).

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

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

Рассмотрим подробнее алгоритм блочного метода. В общем случае четырехугольная зона представляет собой четырехугольник, стороны которого могут быть кривыми второго порядка (рис. 4.5). Как известно, подобную фигуру можно получить как отображение стандартного четырехугольного элемента – квадрата. Каждой точке (x, y) четырехугольника ставится в соответствие точка (, ) квадрата, и наоборот. Локальные координаты (, ) в пределах квадрата (а соответственно и четырехугольника) изменяются от –1 до 1, и, зная, на сколько частей необходимо разделить стороны зоны, легко определить локальные координаты узлов сетки на линиях и внутри зоны (рис. 4.5, б).

Рис. 4.5. Дискретизация четырехугольной зоны блочным методом.

Затем по этим локальным координатам легко найти декартовы координаты узлов сетки внутри исходной зоны:

,

где X_m и Y_m – декартовы координаты вспомогательных узлов (вершин) четырехугольной зоны, {N_m, m=1,…,8} – функции формы для четырехугольного изопараметрического элемента 2-го порядка.

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

Рис. 4.6. Пример разбиения зон блочным методом.

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

Рис. 4.7. Триангуляция методом Делоне. Заштрихованы элементы, которые должны быть удалены (их описанные окружности содержат внутри новый узел).

Различные варианты метода Делоне, как правило, отличаются друг от друга алгоритмами вычисления координат новых внутренних узлов. В данном препроцессоре реализованы два основных алгоритма добавления новых узлов. Первый предусматривает построение “центрального” узла, координаты которого вычисляются с учетом плотности разбиения границ зоны. Затем строятся отрезки, соединяющие центральный узел с углами зоны. Эти отрезки разбиваются на одинаковое число частей, получившиеся при этом узлы соединяются отрезками так, что образуется ряд вложенных фигур, являющихся (приближенно) образами внешней границы при гомотетии относительно центрального узла. Стороны этих вложенных фигур также делятся узлами на части, так, чтобы их длина не отличалась сильно от длины ранее созданных ребер. Преимущество такого алгоритма в его высокой скорости и простоте, а недостатками являются невозможность разбиения сильно вытянутых зон, зон, где необходимо добиться очень сильного изменения плотности сетки, а также зон, форма которых сильно отличается от выпуклой (рис. 4.8).

Рис. 4.8. Этапы дискретизации зоны методом Делоне. а) начальное разбиение (2 конечных элемента); б) добавлен 1-й из узлов будущей сети; в) добавлены все узлы сетки; г) удалены лишние элементы (расположенные вне зоны) и проведена оптимизация расположения внутренних узлов.

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

Все полученные таким образом узлы вносятся в конечно-элементную сетку по алгоритму Делоне, рассмотренному выше. Затем элементы, которые расположены за пределами разбиваемой зоны, удаляются.

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

Рис. 4.9. Локальная оптимизация сетки путем перемещения общего для нескольких элементов узла в барицентр образованной ими фигуры.

Блок проверки качества сетки. На точность конечно-эле­ментного решения в первую очередь оказывает влияние ошибка дискретизации. Ее можно уменьшить, с одной стороны, уменьшая размеры элементов, т.е. увеличивая их число на единицу длины, с другой стороны, увеличивая порядок аппроксимирующих полиномов на каждом конечном элементе. Кроме этой, вполне естественной ошибки, на точность может существенно повлиять форма конечного элемента. Слишком вытянутые, “при­плюснутые” элементы (рис. 4.10) ухудшают обусловленность глобальной системы линейных алгебраических уравнений, что приводит к потере устойчивости решения и, как следствие, к увеличению его погрешности. Вот эту, в общем-то, необязательную ошибку необходимо сводить к минимуму, что обеспечивается получением элементов, форма которых удовлетворяет критерию соразмерности всех сторон и внутренних углов.

Рис. 4.10. “Хороший” (а) и “плохие” (б, в) треугольные конечные элементы.

Критерием качества треугольного конечного элемента является отношение радиусов вписанной и описанной окружностей: k = r/R. Максимально возможным значением этого отношения, как известно, является k = 1/2 (правильный треугольник). Если для рассматриваемого элемента эта величина оказывается меньше 0,15–0,20, то данный элемент признается плохим, и конечно-элементная сетка в соответствующей зоне должна быть удалена и сгенерирована заново.

Препроцессор содержит средства, позволяющие оценить качество сформированной сетки. В случае треугольных элементов, например, на экран может быть выведена цветовая карта, на которой каждый конечный элемент закрашивается определенным цветом в зависимости от величины k. Также пользователь может задать какое-либо предельное значение k₀ из диапазона (0, 1/2), и при этом будут закрашены только те элементы, для которых k < k₀.

Последним этапом работы препроцессора перед запуском решателя является сохранение данных о задаче в виде файлов на диске. Вид уравнения, граничные и другие условия, заданные на линиях (поверхностях), неоднородности среды, геометрия области, другая служебная информация сохраняется в файле taskinfo.xxx, где xxx – расширение имени файла – число от 1 до 999 – идентифицирует задачу. Все параметры конечно-элементной сетки сохраняются в файлах:

xr.xxx, yr.xxx – координаты узлов;

nn1.xxx – списки узлов, образующие конечные элементы;

potent.xxx – привязанные к сетке граничные условия.