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



, О

-

Новосибирский государственный технический университет

г4

На правах рукописи

Рояк Михаил Эммануилович

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

Специальность 05.13.16. - Применение вычислительной техники, математического моделирования и математических методов':.! научных исследованиях (в области технических наук)

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

Новосибирск - 1997

Работа выполнена в Новосибирском государственном техническом университете

Научные руководители : к.т.н. доц. Ю.Г.Соловейчик

к.т.н. доц. Э.П.Шурина

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

научный сотрудник В.Д. Лисейкин

кандидат физико-математических наук И.И.Шабалин

Ведущая организация: ч Вычислительный центр СО РАН

Защита состоится 1997 года в часов на заседании

диссертационного совета Д063.34.03 при Новосибирском Государственном Техническом Университете (630092, Новосибирск-92, пр.К.Маркса, 20).

С диссертацией можно ознакомиться в читальном зале библиотеки НГТУ.

Автореферат разослан 1997 г.

!

• Ученый секретарь /

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

совета (Ук.т.н, доц. Г-П.Чикильдин

Общая характеристика работы.

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

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

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

В проблеме построения конечноэлементиых сеток наиболее важными являются следующие аспекты:

• описание сложной геометрии;

• описание сгущений и разрежений узлов сетки в местах резких изменений градиента решения;

• автоматизация построения сетки по заданному описанию.

Для двумерного случая в научных публикациях можно встретить достаточно много различных подходов к решению проблемы построения конечноэлементиых сеток. Эти подходы обсуждались в работах Jin П., Joe В., Mavriplis D.J., Rebay S., Кузнецова А.Ю., Шайдурова B.B. и др. Последние исследования в этой области направлены в основном либо на повышение эффективности широко известных алгоритмов, например, алгоритма Ватсона построения триангуляции Делоне, либо на автоматизированное построение адаптивных сеток. Заметим однако, что для построения эффективной адаптивной сетки всо- равно необходимо достаточно точно описывать сложную геометрию и строить достаточно близкую к оптимальной первичную сетку, поэтому подходы, использующие адаптивные треугольные сетки, можно рассматривать как еше одну область применения рассматриваемых в данной работе алгоритмов построения эффективных первичных триангуляции.

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

1) методы, основанные на построении трехмерной триангуляции Делоне по произвольному набору узлов;

2) методы, основанные на фронтальном распространении сегки с одновременной генерацией узлов..

Методы первого типа чаще всего применяются в сочетании с описанием области методом конструктивной геометрии. Большинство из них основано на алгоритме, который предложил Watson D.F., а затем усовершенствовали для трехмерного случая Cavendish J.С., Field D.A. и FreyW.H. Очевидным достоинством такого подхода является возможность построения трехмерной триангуляции при практически произвольном набросе узлов (следует отметить, что произвольная генерация узлов в трехмерном случае с заданием необходимых сгущений и разрежений сетки является отдельной, достаточно сложной проблемой, в том числе и из-за трудности визуализации получаемой сетки). Однако триан-

гуляция Делоне в трехмерном случае не только не учитывает внутренние границы области, но и может содержать вырожденные или почти вырожденные тетраэдры, и поэтому такая тетраэдральная сетка не может быть использована в качестве конечноэлементной без ее предварительной коррекции. С учетом затрат на коррекцию сетки время работы алгоритмов построения трехмерной триангуляции Делоне с ростом числа узлов N растет пропорционально Л' \ Таким образом, трехмерная триангуляция Делоне оказывается совсем не оптимальной как с точки зрения качества конечноэлементной аппроксимации, так и с точки зрения затрат вычислительных ресурсов на ее построение при необходимости сколь-нибудь значительного увеличения числа узлов в конечноэлементной сетке.

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

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

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

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

тромагиитных или тепловых нолей любого уровня сложности. Эти алгоритмы реализованы в программном комплексе ТГ;1.МЛ, их эффективность подтверждена решением значительного числа практических задач.

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

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

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

Научная иовнша работы состоит в следующем:

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

2. Предложен метод ускорения алгоритма Ватсоиа построения двумерной триангуляции Делоне, требующий меньшего по сравнению с известными методами объема дополнительной оперативной памяти.

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

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

5. На основе решения большого числа модельных и практических задач выработаны рекомендации по построению близких к оптимальным конечноэлементных сеток.

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

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

Личный вклад. Все результаты, изложенные в диссертации без ссылок на чужие работы, принадлежат лично автору.

Апробация работы. Основные результаты работы были представлены и докладывались на: научном семинаре в техническом университете г.Хемниц (Германия, 1990); Всесоюзной конференции по методам численного решения многомерных нестационарных задач математической физики (Арзамас-16, 1991), Международной геофизической конференции и выставки по разведочной геофизике SEG-EAGO (Москва, 1993); Международном совещании-семинаре по механике реагирующих сред и экологии (Томск, 1994); First Asian Computational Fluid Dinamics Conference (Hotig Kong, 1995); Particle Accelerator Conference and International Conference on High-Energy Accelerators (Dallas,1995); Международной конференции PaCT-95 (Санкт-Петербург, 1995); Международной геофизической конференции и выставке Санкт-Петербург'95; Международной конференции "неклассическая геоэлектрика" (Саратов, 1995); Международной конференции "Математические модели и численные методы механики сплошных сред"' (Новосибирск, 1996); Научно-техническом совещании "Геофизические методы при разведке недр и экологических исследованиях" (Томск, 1996); Международной геофизической конференции "Электромагнитные исследования с контролируемыми источниками" (С.-Петербург, 1996); научных семинарах ВЦ СО РАН, ИВТ СО РАН (Новосибирск). Результаты проведенный научных исследований включались в отчеты по НИР НЭТИ, заключительные отчеты СНИ-ИГГиМСа.

Публикации. По результатам выполненных исследований опубликовано 20 печатных работ.

Структура работы. Диссертация состоит из введения, четырех глав, заключения и списка использованных источников (124 наименования). Работа изложена на 170 страницах, включая 32 иллюстрации.

Основное содержание работы

Глава 1. Описание двумерных расчетных областей н сеток в них

В первой главе рассматриваются методы описания двумерных расчетных областей и средства управления генерацией узлов в них. В п. 1.1 и п. 1.2 описываются Способы и средства ввода и хранения геометрии расчетной области через задание границ ее подобластей.

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

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

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

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

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

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

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

' = 2,п (2)

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

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

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

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

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

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

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

В п. 1.4 предлагается алгоритм автоматического наброса узлов в части расчетной области, ограниченной многоугольником. Алгоритм основан па отображении многоугольника на некоторый четырехугольник, на каждой из сторон которого задано различное число узлов, и последующем переносе сгенерированных на четырехугольнике узлов обратно на многоугольник. При этом координаты узлов I''' =(/','',/''') внутри единичного квадрата для случая, когда на

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

Г>:

I

т + Г

ш+ I

+-1

. У = I.

т+ 1

1т, (3) i

где т - число внутренних узлов на сторожах >=() и ; -1; н0, »,„,1 - число внутренних узлов на сторонах л=0 и л= 1; |р) - целая часть числа р.

Координаты I''1 = , | соответствующих узлов выпуклого четырехугольника рассчитываются по формулам.

Р" = АР?(\ - /")+- /',")(! - г)+г(1 -- Г)К' + ¡К''?.

(4)

где А = й = С = Л = - координаты вершин

четырехугольника, на сторонах АН и С1) которого задано по т узлов, на стороне НС- Пи узлов, на стороне А/) »,„+1 узлов.

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

Пусть на стороне>=0 задано >|/ц внутренних узлов, на сторонезадано*!;, узлов, на стороне л=0 задано Е,и внутренних узлов и на стороне -г=1 задано £| внутренних узлов. Пусть, для определенности, <4,, Ч'(, < 11 выполняется

7

V

либо условие: < , либо, при = , В,. £ Основная идея алгоритма

V, Ч»,

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

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

Пусть Х{ - координата линии, при Пересечении которой количество узлов должно измениться с / на (/+1). Тогда номер к» узла на стороне у=0, который должен использоваться для построения трапеции, можно определить по формуле

Аналогично определяется номер узла кна стороне >--1:

К =[-х|(ч/! +))+оз|. (6)

При этом значением можно вычислить как

+ 1 . " (7)

Однако, если для построения вершин отсекаемой трапеции использовать только формулы (5}-(7), можно не добиться желаемого результата, поскольку отсутствует запрет на получение треугольников. Кроме того, иногда возможно получение сетки с резким изменением шага. Например, если на одном основании трапеции окажется один внутренний узел, а на другом - ни одного, то шаг изменится вдвое. Поэтому необходимо произвести корректировку вычисленных по формулам (5)-(6) значений к\ и к\, учитывая значения к^1 и А:,1'1 для предыдущего отсечения. Эга корректировка дс, может быть выполнена по формуле:

Х- ~ шах( ——— ,х |. (8)

+ ' )

Выполним следующие корректировки:

1 .Если кц < положим к'а = , изменим дг„ используя (8), и пересчитаем к\ по формуле (6).

2.Если к\ < к\л, положим к,' = А,'"1.

3.Если = +1 и к\ = к\А + 2, положим к^ - кЦ' + 2, изменим используя (8), и пересчитаем к[ по формуле (6).

4.Если к'0 = + I и к{ > А,"' + 3, положим к'0 = к'„'' + 2.

5.Если после всех корректировок к'а = А,'/1 или к\ = к*', то текущая трапеция не строится (поскольку она вырожденная).

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

/ = + + /). ^ (9)

Для того, чтобы избавиться от резких переходов шага, можно использовать для генерации узлов внутри трапеции не полученное по формуле (9) число узлов / на боковой стороне, а число яг\ определяемое как:

Ш1пС

+ 1, / • (Ю)

При этом на самой боковой стороне должно быть сгенерировано тт(/и'+,1,0 внутренних узлов.

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

Пусть задан многоугольник, для четырех вершин которого Л =

й = ,Ву|, = Г) = (А<> } задано соответствие вершинам четырех-

угольника, а для вершин четырехугольника - соответствие вершинам единичного квадрата. Будем считать, что на каждой стороне этого единичного квадрата задано число внутренних узлов, равное количеству вершин исходного многоугольника между двумя его вершинами, соответствующими вершинам рассматриваемой стороны единичного квадрата. Для генерации узлов внутри многоугольника можно построить отображение /•', переводящее узлы, сгенерированные в единичном квадрате, на внутренность многоугольника. Будем строить отображение !•' так, чтобы прообразом вершины А многоугольника была точка (0,0), то есть = (0,0); аналогично для остальных вершин многоугольника:

г>(в)"(о.О. '-"'(¿'Ми), гр) = 0,о).

Для описания отображения потребуется параметрическое представление ломаной - образа стороны единичного квадрата. Введем вектор-функцию /'(.г'',*'>•••>■*" определяющую ломаную через параметр

reI0.ll:

.....« (И)

Обозначим через /\ц(') параметрическое представление ломаной с концами А и Н и внутренними узлами х\...,х" (Аналогично введем обозначения 1'ж(1), 1'сп{1), /»,('), /»Л(') и т.д. (при этом, очевидно, выполняется Рдп(/) = ЯВЛ(1 - <) и т.д.). Обозначим также через \АИ\ евклидово расстояние между точками А и В (АВ, = В- А;), и А В - длину ло-

п-1

маной с концами А и В (АВ -\Ах'; + У 1лг|'дг''>'; + \х"В., где - внутренние точки

ломаной). Пусть § = Д.), Л', е[0,1], е{0,1] - точка единичного квадрата. Тогда ее образ = Д} можно рассчитать следующим способом. Вычислим вначале образы проекций .V на стороны квадрата:

йн = Щ =

(12)

и коэффициенты

С _ ^'лВ С _ , _ г _ ^Г" П^Л

5д»- ш пс ~ /Ю •

Рассмотрим следующие ситуации. 1. Две какие-либо противоположные стороны-ломаные являются на самом деле прямыми линиями.

Пусть, для определенности, I АВ\ ~ А В и - СО. Тогда

Л;' ~ (I - + а5,'с, 04)

где

2. Никакие две противоположные стороны-ломаные не являются прямыми линиями.

Обозначим через ОН прямую, проходящую через точки О и Н. Построим точку (?={(/,,{/,), являющуюся пересечением прямых Л".,, и А','^. Л',',». Если прямые не пересекаются (что возможно только в случае сильно невыпуклого многоугольника), то образа точки Л' не существует, и ее необходимо исключить из рассмотрения. В остальных случаях построим еще две вспомогательные точки

К' = (1-«>?,',о + а?^ ,

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

где

Тогда, обозначив ш, = 1НГ>\ и <а,=//Ка|, можно записать:

{/.при ш, = 0 или со, =0

">; Г' ш,+(0г ш,+сог

Рис. 2. Генерация узлов в многоугольнике, близком к четырехугольнику

(20)

На рис. 1-3 приводятся примеры использования разработанных методик и алгоритмов генерации узлов для наброса узлов на некоторые макроэлементы.

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

Рис. 3. Генерация узлов в многоугольнике, блнжом к треугольнику

ку, а также алгоритм выделения всех макроэлементов. Далее обсуждается проблема автоматического построения соответствия многоугольника четырехугольнику или треугольнику для наброса узлов. В п.1.6 описываются различные аспекты' реализации рассмотренных методов в программном комплексе ТЕЬМА.

Глава 2. Построение двумерной конечноэлементной сетки

Во второй главе рассматриваются алгоритмы построения двумерной триангуляции по заданному набору узлов. В п.2.1 анализируется алгоритм Ватсона построения триангуляции Делоне. В п.2.2 предлагается алгоритм коррекции триангуляции Делоне, необходимый для учета внутренних и внешних границ расчетной области. В п.2,3 и п.2.4 рассматриваются методы ускорения алгоритма Ватсона за счет ускорения поиска треугольников и предварительного упорядочивания узлов. В п.2.5 обсуждаются варианты построения первичной триангуляции для инициализации алгоритма Ватсона. В п.2.6 проводится сравнительный анализ скорости построения сетки предлагаемых вариантов реализации алгоритма Ватсона как на регулярных сетках, так и на нерегулярных сетках, использованных для решения практических задач. ^,

В таблице 1 приведены т / . п, . .„

- Таблица 1. Обоощенные данные о вычислительных затра-

обобщенные данные о алия- тах на построение регулярных сеток.

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

Из таблицы I видно, что время работы традиционного алгоритма Ватсона при построении регулярных сеток растет пропорционально к1, а время работы алгоритма с использованием рассмотренных ускорений - лишь немного быстрее, чем теоретически н?улуч-шаемый порог 0(И И).

В таблице 2 приведены данные о вычислительных затратах на построение сетки, содержащей 4600 узлов. Эта сетка изображена на рис.4.

Данные, приведенные в таблице 2, и аналогичные данные, полученные для других нерегулярных сеток, показывают,

Время получения триангуляции (сек)

Число узлов в сетке Без ускорений С ускорениями

Всего Упорядочивание

1024 8 6 1

2500 41 20 4

4900 142 45 7

6400 238 71 14

10000 605 133 26

18225 2044 269 45

Таблица 2. Вычислительные затраты на построение нерегулярной сетки, содержащей 4600 узлов.

Ускорение поиска треугольников - + - +

Время упорядочивания 0 0 17 17

<еек)

Время работы 120 53 К® 40

алгоритма (сек)

Общее время получения трнангулячии (сек) 120 53 136 57

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

Рис. 4. Фрагмент триангуляции расчетной облает для моделирования электродвигателя; эта триангуляция содержит 4600 узлов и 9170 треугольников

Глава 3. Генерация трехмерных конечноэлементных сеток

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

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

образом.

Номер Л'ц рассматриваемого узла сетки определяется по формуле

= (21)

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

(22)

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

Пусть I' = - глобальные координаты узла одного из основных

сечений, а О = (с/,- координаты соответствующего ему узла второго

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

узлам /' и у основных сечений, определяются следующим образом.

Если между двумя основными сечениями задана генерация промежуточ-

ных сечений по прямой, то и + 11 '

Д =

Если же между двумя основными сечениями задана генерация промежуточных сечений по дуге, ось которой определена точками А и В, то

и+1 п+1

(24)

где Й, - точка, которая получилась бы, если бы локальные координаты узда Р были равны заданным локальным координатам узла (), а Л, - точка, которая получилась бы, если бы локальные координаты узла (7 были равны заданным локальным координатам узла Р. Эти точки можно вычислить следующим образом:

Я, = соб(со'<р|е-Щ+ вЦса'<р) [ё-^]*^} = соз((ю1 - - Рj + яп({и>' -1 )Щр

(25)

(26)

где .V = А) - направление оси поворота, (5 и Р - проекции на ось по-

ворота точек У и Р соответственно^ которые можно вычислить по формулам ¿М + р-Л],^, (27)

= Л + (28)

Ф - угол между точками {? и Р, определяемый как

Ф = агс1ц

( Р - О - * ^

(29)

а а' -

вес /-ого узла, вычисляемый, по формуле

(и-/ + !)(»+!), с=1

(о1 = (30)

(с'"1 -с"+')/(1-с"+'), с*1 .

В п.3.5 описывается реализация рассмотренного подхода к построению трехмерных конечноэлементных сеток в комплексе программ ТЕЬМА.

Глава 4. Построение сеток при решении некоторых модельных и практических задач методом конечных элементов

В п.4.1 обсуждаются особенности построения конечноэлементных сеток при моделировании стационарного электрического поля для задач электроразведки в проводящей среде с учетом обсадной колонны труб. В п.4.1.1 приводится пример осесимметричной задачи, для которой триангуляция Делоне не является оптимальной, предлагается способ ее коррекции, основанный на описанном в п.2.2 диссертационной работы алгоритме и сравниваются результаты расчета дня одного и того же набора узлов при использовании триангуляции Делоне и откорректированной триангуляции. В п.4.1.2 приводится пример влияния распределения узлов в трехмерной сетке на скорость сходимости итерационного процесса решения конечноэлементной СЛАУ, аппроксимирующей трехмерную краевую задачу. В п.4.1.3 приводится пример построения конечно-элементной сегки для расчета стационарного электрического поля в проводящей среде при наличии нескольких скважин с обсадными колоннами труб.

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

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

- 1 дЛ + оГ ^ +егаап = .7 , (31)

Мо V з/ )

-а|у(авгааг)-^у|а^| = -й\'7 , (32)

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

Для задачи расчета осесиммегричной части поля в диссертационной работе сформулированы рекомендации по построению конечноэлементных сеток, приводится пример расчета изменения осесимметричного электромагнитного поля во времени. Здесь же рассматривается пример использования предлагаемых в третьей главе методов построения трехмерной сетки для расчета поля влияния трехмерного объекта. Сечения этой сетки плоскостями г=0 и у=0 приводятся на рис.5-6.

7 ммж

■жш

Рис. 5. Сечение плоскостью г=0 теграэдраль- Рис. 6. Сечение плоскостью>"=0 тетраэдральной сетки, использованной для расчета поля ной сетки, использованной для расчета поля влияния трехмерного объекта. влияния трехмерного обьекта (масштабы по

осям х и г разные).

В п.4.4. описываются специальные средства построения сеток, созданные для моделирования нестационарных электромагнитных полей в технических конструкциях. Для вычисления этих полей используется вариационная постановка, предложенная Соловейчиком Ю.Г. Эта постановка требует для своей корректности задания в расчетной области специальных поверхностей-перемычек, прослоек-, перемычек и токовых сечений. В диссертационной работе приводится пример построения сетки в одной из таких конструкций. Общий вид расчетной области для моделирования электромагнитного поля в этой конструкции приведен на рис.7. С'сченне сетки плоскостью х=0 приводится на рис.8. Утолщенными линиями выделены внутренние границы расчетной области и заданная поверхность-перемычка. Зэшгрнховано поперечное сечение обмотки, по которому вычислялся полный ток (эта поверхность одновременно является частью поверхности симметрии д~0).

|

мщ

шщ

Ж^УМ

Рис. 7. Общий вид расчетной области для моделирования нестационарного электромагнитного поля в технической конструкции.

/ \ \ 1 V \ / /' / \ / / \ / /

/ ,/ \ \ '/ \ \ \ I / \ / / /

\ / \ / \ / X

\ / \ \ / \ / \ \ [1 /; 1 / / / У

/ / \ / / N \ \ \ / / \ / \ 1 \ ТА- . ш 1 / / / / /

/ \ \ / / / / /

V / \ \ / / \ 1 ж / / / /

\ \ у -V __-—"

/ / \ / \ \ 7 \ / / / \ \ /

\ / \ \ У \ У / \ / //

/ Ч \ / \ 7 \ \ / \ \/ / \ /

/ \ / \ \ У \ / \ / \

Рис. 8. Сечение конеч ^элементной сетки плоскостью г=0.

Основные результаты работы

На основе проведенных в диссертационной работе исследований получены следующие теоретические и практические результаты.

1. Разработаны методы описания двумерной расчетной области, позволяющие задавать геометрию практически произвольной сложности и эффективно управлять размещением узлов. Созданы алгоритмы автоматического и полуавтоматического наброса узлоз по заданному описанию: Предложенные методы и алгоритмы реализованы в двумерном препроцессоре пакета ТЕЬМА и использовались при решении значительного числа практических задач.

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

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

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

4. На основе метода тиражируемых сечений разработаны средства описания тетраэдральной сетки в сложной трехмерной геометрии, позволяющие эффективно управлять перемещением узлов. Предложен способ формирования сложных поверхностей с использованием вспомогательной триангуляции - маски. Разработан однопроходный алгоритм построения согласованной конечно-элементной сетки по заданным сечениям, не использующий никакой дополнительной информации о триангуляции поверхности, кроме списка узлов и списка треугольников. Разработанные средства и алгоритмы реализованы в трехмерном препроцессоре пакета ТЕЬМА.

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

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

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

Основное содержание диссертации отражено в следующих работах

1. Шурина Э.П., Соловейчик Ю.Г., Соломахин В.И , Рояк М.Э. Триангуляция сложных областей в задачах моделирования теплового состояния электрических машин. П Машинные методы оптимизации, моделирования и планирования эксперимента. - Новосибирск, НЭТИ, 1988. - с.86-88.

2. Создание программного комплекса моделирования теплового состояния криостатнрующего устройства : Отчет о НИР/ НЭТИ; руководитель Э.П.Шурина . - № 01.9.10 019648; Инв. № 02.9.10 026252. - 1991. Hobqch-бирск. - 54с.

3. Шурина Э.П., Соловейчик Ю.Г., Рояк М.Э. Особенности моделирования нелинейных физических процессов в трехмерных областях. // Вопросы атомной науки и техники. Серия Математическое моделирование физических процессов. - М., 1992, Вып.З. - с.86-87

4. Шурина Э.П., Соловейчик Ю.Г., Рояк М.Э. Моделирование электромагнитных полей в трехмерных областях. // Вычислительные технологии. - Новосибирск, Институт вычислительных технологий СО РАН, 1993, Т.2, №6 - с.48-53.

5. Моисеев B.C., Соловейчик Ю.Г., Рояк М.Э. Математическое моделирование сложнопостроенных сред. II Сборник рефератов №2 Международной геофизической конференции и выставки по разведочной геофизике SEG-EAGO. -М., 1993.-с. 15.

6. Шурина Э.П., Соловейчик Ю.Г., Рояк М.Э. Математическое моделирование физических полей, обусловленных локальными возмущениями. // Вычислительные технологии. - Новосибирск, Институт вычислительных технологий СО РАН, 1994, Т.З, №8. - с.143-147.

7. Шурина Э.П., Соловейчик 10.Г., Рояк М.Э. Моделирование физических полей в трехмерных объектах. // Сопряженные задачи физической механики и экология: Тезисы докладов международной конференции.-Томск, 1994.—1с.

8. Моисеев B.C. Соловейчик Ю Г. Рояк М.Э. Тригубович Г.М. Влияние диэлектрической проницаемости среды на распространение электромагнитной волны. // НЕКЛАССИЧЕСКАЯ ГЕОЭЛЕКТРИКА: Материалы международной конференции (28 августа - 1 сентября 1995г., г. Саратов. - Саратов: Нижневолжскнй НИИ геологии и геофизики, 1995. - стр.60-61.

9. Моисеев B.C., Соловейчик Ю.Г.; Рояк М.Э., Тригубович Г.М. Математическое моделирование электромагнитных полей в сложных средах. // Международная геофизическая конференция и выставка Санкт-Петербург'95 (10-13 июля 1995г): Тезисы докладов. - С.-Петербург, 1995г. - Т.2, №18.4

10. Grudiev A., Royak М., Shurina E.,.Soloveichik Yu, TiunovM., Vobly P.. Mathematical Simulation of 3-D Magnitostatic Fields Using the Complex of Programs MASTAC. // Abstracts of AMCA?95. - Novosibirsk: NCC Publisher, 1995. -P. 131 -132

1 l.Rojak M., Shurina E., Soloveichik Yu. and Malyshkin V. Parallelization of Computer Code MASTAC Three-Dimensional Finite Elements Method Implementing // Proceedings of PaCT-95 Lecture Notes in Computer Science. - Germany: Springer, 1995,-P.304-313.

12. Shurina E.P., Soloveichik Yu.G., Royak M.E., Vobly P.D., Grudiev A.V., Tiunov M.A. Three-Dimensional Non-linear Magnitostatic Fields Mathematical Modelling // The Scientific Conference on Use of Research Conversion Results in the Sibe-

rian Institutions of Higher Education for International Cooperation (SIBCONVERS-95). - Tomsk, 1995. - P.30.

13.Shurina E.P., Soloveitchik J.O., Royak M.E. Three-dimensional Fields Modeling on Irregular Mesh Using Finite Elements Method // Proceedings of the First Asian Computational Fluid Dynamics Conference 16-19 January. - V.3, Hong Kong, 1995. - P.1225-1226

14. Моисеев B.C., Соловейчик Ю.Г., Рояк М.Э., Тригубович Г.М. Новое в математическом моделировании электромагнитных полей искусственных источников тока. // Международная геофизическая конференция "Электромагнитные исследования с контролируемыми источниками": тезисы докладов - С.-Петербург, 1996. - с.22-23.

15. Моисеев B.C., Соловейчик Ю.Г., Рояк М.Э., Васильев А.В. Математическое обеспечение задач электроразведки для сложнопостроенных сред. // Научно-техническое совещание "Геофизические методы при разведке иедр и экологических исследованиях": сборник материалов (докладов).-Томск, 1996 - с.65-66.

16. Рояк М.Э., Соловейчик Ю.Г. Алгоритмы построения нерегулярных тре-. угольных и тетраэдральных сеток. // Сб. научных трудов НГТУ, - Новосибирск, НГТУ, 1996г., №2(4). - с.39-46.

17. Соловейчик Ю.Г., Рояк М.Э. Математическое моделирование трехмерных нестационарных электромагнитных полей в электротехнических устройствах // Труды третьей международной научно-технической конференции "Актуальные проблемы электронного приборостроения" АПЭП-96 в 11 томах./Т.6 Моделирование и вычислительная техника. - Новосибирск, НГТУ, 1996г.-с.66-69

18. Шурина Э.П., Соловейчик Ю.Г., Рояк М.Э. Вычислительные схемы решения трехмерных нелинейных магнитостатических задач. // Тезисы докладов международной конференций "Математические модели и численные методы механики сплошных сред". - Новосибирск, 1996 - с.529.

19. Шурина Э,П., Соловейчик Ю.Г., Рояк М.Э. Решение трехмерных нелинейных магнитостатических задач с использованием двух потенциалов. Новосибирск. 1996. - 28 с. (Препринт / РАН. Сиб. отд-ние. ВЦ ; № 1070).

20. Rojak М., Shurina Е., Soloveichik Yu. and Grudiev A., Tiunov M., Vobly P. MAS-TAC - new code for solving three-dimensional non-linear magnetostatic problems // Proceedings of Par-95 (Particle Accelerator Conference and International Conference of High-Energy Accelerators).-Vol.4, USA, Texas, Dallas: ШЕЕ, 1996.-P.2291 -2293. ■

21. Soloveichik Y.G., Royak M.E. The computing schemes of non-stationary electromagnetic fields FEM modeling in mediums with three-dimensional inhomogeneity. // Bulletin of the Novosibirsk Computing Center. Series: Computer Science. - Novosibirsk: NCC Publisher, Issue 4, 1996. - P.55-78

22. Разработка технологии площадных работ с многоканальной аппаратурой при поисках залежей углеводородов в Сибири. Заключительный отчет. Отв. исполнитель Моисеев B.C. - СПИИГГиМС, Новосибирск, 1996г. - 145с.