Содержание к диссертации
Введение
ГЛАВА 1. Подготовка данных для начально-краевых задач 18
1.1. Применение САПР для подготовки данных начально краевых задач 18
1.1.1. Основные модели представления форм 19
1.1.2. Обмен данными между системами геометрического моделирования 26
Формат IGES 27
Формат DXF 30
Формат STEP 31
1.2. Модель NURBS 39
1.2.1. Свойства В-сплайн базиса 40
1.2.2. В-сплайн кривые 41
1.2.3. В-сплайн поверхности 47
1.2.4. В-сплайн объемы 49
1.2.5. Нахождение прообраза точки В-сплайн кривой 54
1.2.6. Нахождение прообраза точки В-сплайн поверхности 57
1.3. Геометрически комплексы и импорт данных 60
1.3.1. Введение в геометрически комплексы 60
1.3.2. Преобразование данных при импорте в геометрические комплексы 64
1.4. Алгоритмы генерации расчётных сеток 66
1.4.1. Алгоритм построения двумерных гибридных сеток 67
1.4.2. Алгоритм триангуляции составной В-сплайн поверхности 71
1.5. Краткие выводы 75
ГЛАВА 2. Разработка методики расчёта притока флюида к горизонтальной скважине 77
2.1. Основные положения метода конечных элементов 77
2.1.1. Методы Ритца и Галеркина как основа МКЭ 79
2.1.2. Классы элементов, используемые в программном комплексе 80
2.1.3. КЭ аппроксимация уравнения Пуассона 81
2.1.4. КЭ аппроксимация уравнения диффузии 86
2.1.5. Подготовка расчётной NURBS геометрии 86
2.2. Моделирование пластовых систем и скважин 89
2.2.1. Закон фильтрации 91
2.2.2. Задача однофазной фильтрации 94
2.2.3. Интерпретация граничных условий в задачах фильтрации 96
2.2.4. Эффективный радиус скважины и проблема разномасштабности 99
2.2.5. Алгоритм макроблока 101
2.2.6. Аналитические решения и модель начального приближения .110
2.3. Краткие выводы 116
ГЛАВА 3. Модельные расчёты фильтрации и примеры работы основных алгоритмов 118
3.1. Тестовые модельные примеры работы геометрических алгоритмов 118
3.1.1. Нахождение прообраза NURBS кривой 118
3.1.2. Нахождение прообраза NURBS поверхности 120
3.2. Модельные примеры расчета задач фильтрации 123
3.2.1. Решение Лейбензона. Анализ характера течения вблизи скважины 123
3.2.2. Плоскорадиальное течение 126
3.2.3. Трехмерная однофазная фильтрация на призматической сетке 127
3.2.4. Трехмерная однофазная фильтрация на тетраэдральной сетке 129
3.3. Пример расчёта теплопроводности на NURBS элементах 130
Заключение 132
Литература
- Обмен данными между системами геометрического моделирования
- Нахождение прообраза точки В-сплайн кривой
- Классы элементов, используемые в программном комплексе
- Нахождение прообраза NURBS поверхности
Введение к работе
Актуальность проблемы
Методы и программы для численного решения начально-краевых задач математической физики развиваются на протяжении многих десятилетий. Развитие компьютерных технологий сделало вычислительный эксперимент мощным и эффективным способом решения проблем прикладной математики. Появление новых математических моделей и численных методик, рост объема и сложности современных вычислений делают необходимым создание новых средств автоматизации построения компьютерных моделей, подготовки исходных расчётных данных и решения начально-краевых задач.
Алгоритмы решения начально-краевых задач интегрированы во многие систем автоматизированного проектирования (CAD системы), например, в SolidWorks (), Unigraphics (), T-Flex (), Salome1 (), ИСПА (). Существуют многочисленные коммерческие пакеты прикладных программ, применяемые для исследований и расчётов, такие как Fluent, CFX, StarCD, Cosmos, FlowVision и другие. Всё большее распространение получают программные разработки с открытым кодом, в частности, программные реализации КЭ методов, например, GetFEM (), OOFEM (), FreeFEM++ (), а также КО методов, например, EULER (). В ИПМ им. М.В. Келдыша РАН и ИММ РАН для научных и прикладных исследований были созданы специализированные программные комплексы, такие как
Распространяется в соответствии с GNU LGPL лицензией.
ВВЕДЕНИЕ 6
TEK0H[11], ACTA [12], РАЗРЯД [14], GIMM [16-21], MARPLE [23,24] и другие.
Использование современной вычислительной техники для решения задач математического моделирования делает необходимым разработку и внедрение комплексного подхода к их решению. Программные средства, разрабатываемые для этой цели, должны обеспечивать поддержку различных видов вычислительной работы таких, как построение геометрической модели области решения задачи, генерацию расчетной сетки, составление системы сеточных уравнений соответственно принятой математической модели, выполнение расчёта и обработку полученных данных. К программным средствам автоматизации научных исследований предъявляются требования робастности и устойчивости при изменении входных параметров математических моделей в широком диапазоне значений, а также при переходе от геометрически простых расчетных областей к сложным, многосвязным, разномасштабным областям.
Математические модели физических процессов обычно описываются дифференциальными уравнениями в частных производных, а также интегральными и интегро-дифференциальными уравнениями. Для решения прикладных задач в соответствующих постановках используются, в основном, численные методы, из которых наибольшее распространение получили метод конечных элементов[100-110], конечных разностей [11,13] и граничных элементов [111]. В решении общей вычислительной задачи с использованием вышеуказанных численных методов можно выделить три этапа: предобработку, расчёт и постобработку. Первый этап состоит в задании геометрии расчётной области и определении начальных и граничных условий. Далее производится построение расчётной сетки для заданной геометрии. Выбор типа сетки определяется характером решаемой задачи, геометрией области и требуемым объемом вычислений. Следующий этап - собственно расчёт, результатом которого являются некоторые распределения физических величин по элементам сетки. Заключительный этап предполагает постобработку результатов расчёта, нахождение физических и вспомогательных геометрических характеристик, а также их визуализацию на сетке. Сетка, используемая для визуализации,
ВВЕДЕНИЕ 7
довольно часто не совпадает с той, что используется для расчёта.
Первый и третий этапы являются вспомогательными по отношению к основной задаче решения уравнений математической физики. Однако, в задачах, определенных в пространственных областях сложной формы, вопросы подготовки и анализа геометрических данных настолько усложняются, что превращаются в самостоятельные проблемы, для эффективного решения которых требуются специальные средства автоматизации. Разработка таких средств является одной из проблем, исследуемых в данной работе.
Возможность задания реальной геометрии изучаемых объектов -важнейшее требование к прикладному программному обеспечению. Данную задачу удобно решать с помощью специализированных CAD систем. Одним из распространенных методов описания пространственных тел, задаваемых с помощью параметрических кривых и поверхностей, является метод, основанный на граничном представлении. Граничное представление удобно применять для подготовки данных при решении начально-краевых задач математической физики, т.к. сами границы в этом случае определяются функционально. Выбор вида параметрического представления производится из соображений универсальности. В качестве такого рода геометрических примитивов удобно использовать NURBS кривые и поверхности, позволяющие описывать объекты как с криволинейными, так и с прямолинейными границами. .
NURBS широко применяются для задания кривых и поверхностей во многих современных CAD системах. С помощью В-сплайнов можно представить достаточно сложные поверхности, что позволяет построить гибкую систему моделирования форм [1-4]. Примерами систем, поддерживающих работу с В-сплайнами, являются Rhinoceros (), T-Flex, Unigraphics и др. Модели на основе NURBS используются также в универсальных форматах обмена геометрическими данными, например, они реализованы в известном формате STER Разработка средств импорта топологической и геометрической информации из STEP является одной из проблем, решаемых в данной диссертационной работе. При задании геометрии также определяются
ВВЕДЕНИЕ 8
граничные условия, каждое из которых сопоставляется определенной геометрической сущности. В частности, в двумерном случае условия ставятся в соответствие контурным кривым, а в трехмерном - граничным поверхностям. Сами граничные условия имеют произвольную структуру, что весьма удобно для решения различного рода задач.
Следующий этап подготовки данных для численного моделирования -это генерация сетки, которая производится с использованием информации о геометрии модели, заданной NURBS. В зависимости от характера решаемой задачи с использованием NURBS параметризации можно строить как структурированные, так и гибридные сетки. В рамках граничного представления дискретизация осуществляется иерархически в порядке возрастания топологической размерности элементов модели. Для двумерной расчётной области сначала производится разбиение кривых, дающих систему контуров. Соответственно, один из этапов построения 3D сетки состоит в генерации поверхностной сетки. Дискретизация поверхности может быть выполнена с помощью генерации сетки в пространстве параметров с последующим её отображением в физическое пространство. При этом поверхностная сетка должна удовлетворять определенным критериям качества [6]. Для удовлетворения данных критериев генерация сетки в параметрическом пространстве должна выполняться с учетом кривизны поверхности.
Существенную сложность представляет дискретизация составной поверхности, состоящей из NURBS сегментов, содержащих дыры и сложные границы. Эта задача решается построением триангуляции. Особенность 3D моделей, представляемых CAD системами, состоит в том, что параметризация поверхностей не совпадает с параметризацией лежащих на них кривых. Это приводит к необходимости построения алгоритмов поиска прообразов точек моделируемой NURBS поверхности [22]. Данные алгоритмы реализованы в рамках общего алгоритма дискретизации составной поверхности. Для построения треугольной сетки на поверхности при этом используются алгоритмы двумерной триангуляции в сочетании с параметрическим отображением. Особенность описываемого в данной работе подхода заключается в том, что предлагаемый алгоритм дискретизации
ВВЕДЕНИЕ
позволяет строить качественные треугольные сетки на составных поверхностях, описываемых В-сплайн сегментами сложной формы, и обеспечивает адаптацию сетки к особенностях геометрии поверхности. Методики дискретизации В-сплайн поверхностей также рассмотрены в работах [37,54-59].
Существуют многочисленные программные разработки области анализа геометрии и генерации сетки, в основном осуществляющие построение нерегулярных треугольных и тетраэдральных сеток. Постороение нерегулярных сеток осуществляется полностью автоматически для почти произвольной геометрии, в том числе неодносвязных областей, содержащих дыры. Некоторые из этих систем осуществляют дискретизацию областей заданных ограничениями в виде множества контуров (в 2D) или замкнутым множеством поверхностей (в 3D), т.е. используют в основном собственные форматы импорта данных. Можно отметить несколько таких систем в открытом (free) доступе:
Triangle ( quake/triangle.html),
EasyMesh (),
QHull (),
Tetgen ().
Другие системы осуществляют также импорт геометрии из других форматов таких, как STL, IGES, STEP и имеют более сложные структуры данных, на которых осуществляется процесс дискретизации. К ним можно отнести, например, NetGen () и GMSH ().
В данной работе для генерации объемной тетраэдральной сетки используется метод Делоне с ограничениями, реализованный в системе Tet-Gen. В качестве входных данных для TetGen используется поверхностная сетка, построенная описанным выше алгоритмом триангуляции составной В-сплайн поверхности.
Основным этапом решения начально-краевой задачи является расчёт, осуществляемый после задания геометрии, построения сетки, определения начальных и граничных условий. Алгоритм расчёта основан на методе конечных элементов. Для представления криволинейной геометрии и
ВВЕДЕНИЕ 10
повышения точности вычислений в условиях неоднородности среды применяются элементы высокого порядка. В рамках общего КЭ подхода возможно также использование NURBS элементов. NURBS технология позволяет использовать элементы высокого порядка в сочетании с удобством геометрического описания. Длительность вычислений зависит от объема сетки, типа элементов и постановки задачи. На форму, размеры и взаимное расположение элементов сетки накладываются определенные ограничения, обусловленные требованиями точности и сходимости соответствующих разностных и конечно-элементных схем [110]. Использование для расчёта NURBS геометрии, импортированной из файла формата STEP и обработанной в соответствии с требованиями КЭ анализа, не требует дополнительной дискретизации расчётной области.
В связи с вышесказанным, разработка интегрированной программной системы для численного решения трехмерных начально-краевых задач является особенно актуальной проблемой. Программный комплекс должен включать универсальные средства описания и дискретизации пространственных расчетных областей сложной формы, задания граничных условий различных типов, предоставлять возможность использования конечно-элементных аппроксимаций различных типов, гарантировать точность и надежность получаемых численных решений. При создании современных программных продуктов таких, как GIMM[16-21] и MARPLE [23, 24], используется современный объектно-ориентированный подход, позволяющий снизить затраты на поддержку и модификацию программных модулей.
В последние годы наблюдается резкое увеличение потребления нефтепродуктов и постепенное истощение запасов углеводородов. В связи с этим разрабатываются месторождения, которые ранее считались неперспективными ввиду сложности необходимых для этого технологий, к которым относится прокладка наклонных, горизонтальных и многоствольных скважин. Всё более актуальной становится проблема компьютерного моделирования пластовых течений, в частности, построение моделей притока флюида к скважине для оценки её дебита. Такие модели позволяют на начальном этапе проекта разработки углеводородной
ВВЕДЕНИЕ 11
залежи учесть геометрические и физические характеристики пласта, а также оптимизировать расположение скважин с целью улучшения их эксплуатационных характеристик. Разработанные для моделирования притока флюида программные средства тестировались на задачах расчёта продуктивности нефтяных скважин в условиях, приближенных к реальным условиям добычи углеводородов.
Известно, что характерные размеры пласта на несколько порядков больше размеров скважин. Разномасштабность является важной проблемой, затрудняющей расчёт притока флюида. Для решения данной проблемы в диссертации предложено специальное граничное условие, соответствующее скважине. Контур постановки данного условия при решении уравнений фильтрации отодвигается от границ скважины, благодаря чему отпадает необходимость использования сеток, адаптированных к поверхности скважины, что позволяет понизить вычислительную сложность задачи. Пространство между скважиной и контуром постановки граничного условия называется макроблоком. В области макроблока делаются предположения о характере течения флюида, позволяющие заменить численный расчёт в его пределах аналитическим решением. Тем самым, область решения делится на две части: макроблок и внешнюю часть. Численный расчёт производится только во внешней части, что позволяет существенно снизить количество сеточных узлов по сравнению с методикой расчёта фильтрации в полной области с дискретизацией до границы скважины.
В связи с вышеизложенным, целями диссертационной работы являются Создание интегрированного программного комплекса для прикладных исследований процессов, описываемых в рамках двумерной или трехмерной начально-краевой задачи. Программная система должна обеспечивать полный цикл численного решения нестационарной нелинейной задачи, т.е. выполнять следующие функции: импорт геометрической модели расчетной области в виде совокупности NURBS примитивов из файла формата STEP АР 214 в локальные структуры данных, задание начальных, граничных условий, а также набора выходных параметров, генерацию гибридной сетки конечных
ВВЕДЕНИЕ 12
элементов, формирование системы сеточных уравнений и её решение, обработку результатов расчетов. Разработка и реализация математической модели, решающей проблему разномасштабности размеров скважины и пласта для расчёта продуктивности нефтяных скважин при нелинейном законе фильтрации Форхгеймера. Модель должна учитывать реальную геометрию нефтяного коллектора и скважины.
В соответствии с целью исследования были поставлены следующие конкретные задачи:
Реализовать алгоритмы импорта NURBS геометрии из файлов формата STEP;
Реализовать алгоритмы вычисления В-сплайн кривых, поверхностей и объемов;
Разработать и реализовать алгоритмы поиска прообраза В-сплайн кривых и поверхностей на основе поиска начального приближения с последующей доводкой итерационным методом;
Реализовать алгоритмы построения двумерных, поверхностных и гибридных расчётных сеток на базе NURBS геометрии, импортированной из STEP;
Реализовать построение объемных сеток на основе сгенерированных поверхностных с использованием существующих алгоритмов тетраэдризации.
Разработать и реализовать алгоритм объединения сечений для построения призматических сеток с использованием NURBS параметризации.
Исследовать возможности NURBS для построения блочно-структурированных сеток.
Реализовать метод конечных элементов для решения задач нестационарной теплопроводности и фильтрации.
Исследовать возможность применения NURBS элементов для проведения КЭ расчётов.
ВВЕДЕНИЕ 13
10. Реализовать алгоритм макроблока для решения задачи оценки продуктивности нефтяной скважины при нелинейном законе фильтрации Форхгеймера в двумерных и трехмерных расчётных областях.
Научная новизна:
Разработаны новые алгоритмы и программные средства импорта геометрической и топологической информации в структуры данных, используемые для численного моделирования.
Реализована новая технология подготовки данных и расчёта задач с криволинейными границами на основе использования NURBS элементов, позволяющих повысить точность геометрического описания и управлять порядком аппроксимации. NURBS технология является дальнейшим развитием суперэлементного подхода. Импорт NURBS элементов производится непосредственно из файлов формата STEP.
Разработаны следующие алгоритмы генерации сеток: алгоритм триангуляции поверхностей, состоящих из совокупности NURBS сегментов сложной формы; алгоритм генерации гибридных сеток, объединяющий методы параметрических отображений и построения триангуляции; алгоритмы и программные средства генерации трехмерных сеток, объединяющие кинематические методы, методы параметрических отображений и построение тетраэдризаций.
Реализована математическая модель притока флюида к нефтяной скважине в трехмерной постановке, решающая проблему разномасштабности характерных размеров скважины и нефтяного пласта. На основе построенной модели проведены численные исследования дебита скважины и характера течения в зависимости от физических и геометрических характеристик пласта. В численном эксперименте показана эффективность данной модели.
Научная и практическая ценность В разработанном программном комплексе реализована современная технология подготовки данных и новый метод расчёта с использованием NURBS, что позволяет облегчить
ВВЕДЕНИЕ 14
процесс численного анализа начально-краевых задач. Данный программный комплекс может быть применен для исследования широкого класса проблем естествознания и техники, формулировки которых приводят к смешанным начально-краевым задачам. Модульная структура программной системы позволяет использовать отдельные её блоки в составе других интегрированных систем САПР и АСНИ, а также программных комплексов GIMM и MARPLE. Созданный программный комплекс использовался для исследования трехмерных нестационарных тепловых процессов в нелинейных средах, а также для анализа продуктивности нефтяных скважин. Разработанные программные средства позволяют учитывать реальную геометрию нефтяных скважин и пластов, а также, в случае нелинейного закона фильтрации, выделять области наиболее существенного его влияния в практически важных случаях.
Положения, выносимые на защиту:
Технология подготовки данных на основе NURBS для решения начально-краевых задач.
Комплексы прикладных программ подготовки исходных данных и КЭ расчёта начально-краевых задач для параболических уравнений.
Алгоритм КЭ расчёта на NURBS элементах. Алгоритмы импорта исходной NURBS геометрии и подготовки дискретной модели расчётной области для расчёта методом конечных элементов. Применение NURBS технологии для вычислений по методу суперэлементов.
Модель прискважинного макроблока для моделирования течения флюида (нефти, газа) при нелинейном законе фильтрации Форхгеймера.
Достоверность результатов обеспечена обоснованностью применяемых математических моделей, использованием современных методов решения начально-краевых задач, а также верификацией полученных моделей путем сравнения с аналитическими оценками и результатами других авторов.
ВВЕДЕНИЕ 15
Апробация работы Результаты работы докладывались и обсуждались на следующих конференциях:
Международная конференция "Parallel CFD 2004", Las Palmas de Gran Canaria, Испания, 24-27 мая, 2004 г.
Всероссийская Научно-Техническая Конференция "Параллельные вычисления в задачах математической физики", 21-25 июня, 2004 г., г. Ростов-на-Дону.
Международная конференция "Parallel CFD 2005". University of Maryland, College Park, Maryland, США, 24 - 27 мая, 2005 г.
International Conference on matrix methods and operator equations, Москва, ИВМ PAH, 20-25 июня, 2005 г.
Всероссийская научная конференция "Научный сервис в сети Интернет: многоядерный компьютерный мир",24-29 сентября 2007 г., г. Новороссийск.
Совместный семинар ИММ РАН и кафедры математического моделирования МФТИ под руководством проф. д.ф.-м.н. Е.И. Леванова, 18 октября 2007г.
Основное содержание диссертационной работы изложено в следующих публикациях:
Б. Н. Четверушкин, В. А. Гасилов, С. В. Поляков, М. В. Якобовский, И. В. Абалакин, Е.Л. Карташева, И.В. Попов, П. С. Кринов, С. А. Суков, В. Г. Бобков, А. С. Минкин. Пакет прикладных программ GIMM для решения больших задач гидродинамики на многопроцессорных вычислительных системах. - Сборник трудов Всероссийской научно-технической "Параллельные вычисления в задачах математической физики". Изд-во Ростовского Гос.Университета "ЮГИНФО РГУ", 2004, с. 141-158.
В. N. Chetverushkin,V. A. Gasilov, S. V. Polyakov,M. V. Iakobovski, Е. L. Kartasheva, I. V. Abalakin, I. V. Popov, N. Yu. Romanyukha, S. A. Sukov and A. S. Minkin. CFD Software Project GIMM. Study of Hy-drodynamic Problems Via Parallel Computing. - Parallel Computational Fluid Dynamics. Multidisciplinary Applications/Proc. Of the Parallel CFD
ВВЕДЕНИЕ 16
2004 Conference, Las Palmas de Gran Canaria, Spain (May 24-27, 2004). -Amsterdam - Elsevier B.V. 2005, p. 339-344.
Б. H. Четверушкин, В. А. Гасилов, С. В. Поляков, Е. Л. Карташева, М. В. Якобовский, И.В. Абалакин, В. Г. Бобков, А. С. Болдарев, С. Н. Болдырев, С. В. Дъяченко, П. С. Кринов, А. С. Минкин, И. А. Нестеров, О. Г. Ольховская, И. В. Попов, С. А. Суков. Пакет прикладных программ GIMM для решения задач гидродинамики на многопроцессорных вычислительных системах. - Математическое моделирование, 2005г., том 17, номер б, с. 58-74.
В. N. Chetverushkin, V. A. Gasilov, S. V. Polyakov, М. V. lakobovski, Е. L. Kartasheva, A. S. Boldarev, A. S. Minkin. Data Structures and Mesh Processing in Parallel CFD Project GIMM. - ParCo 2005, International Conference on Parallel Computing, Abstracts, p.41.
B. N. Chetverushkin, V. A. Gasilov, S. V. Polyakov, M. V. lakobovski, E. L. Kartasheva, A. S. Boldarev, A. S. Minkin. Data Structures and Mesh Processing in Parallel CFD Project GIMM. - ParCo 2005: Current and Future Issues of High-End Computing, John von Neumann Institute for Computing, p.351.
B. Chetverushkin, V. Gasilov, M. Yakobovski, S. Polyakov, E. Kartasheva, A. Boldarev, I. Abalakin and A. S. Minkin. Unstructured mesh processing in Parallel CFD project GIMM. - Parallel Computational fluid dynamics: Theory and applications/Proc. of the Int. Conf. Parallel CFD 2005 Conference, College Park, MD, USA /Ed. by A. Deane, A.Eger, J. McDonald, N. Satofuka, G.Brenner, D.R.Emerson, J.Periaux, D. Tromeur-Dervout. Elsevier Publ, Amstedam, Oxford, 2006, p. 501-508.
E. Л. Карташева, А. С. Минкин, В. А. Гасилов Метод триангуляции составных поверхностей, состоящих из В-сплайн сегментов сложной формы. - Математическое моделирование, 2007 г., том 19, номер 10, с. 44-60.
Б. Н. Четверушкин, В. А. Гасилов, С. Н. Болдырев, А. С. Болдарев, Е. Л. Карташева, О. Г. Ольховская, С. В. Дьяченко, А. С. Минкин, С. В. Поляков, М. В. Якобовский. Объектно-ориентированное программирование и организация структур данных для параллельных
ВВЕДЕНИЕ 17
расчётов комплексных задач гидродинамики. - Приложении к журналу "Открытое образование", с. 40-41. 9. В. А. Гасилов, С. Н. Болдырев, А. С. Болдарев, Е. Л. Карташева, О. Г. Ольховская, С. В. Дьяченко, А. С. Минкин. Об использовании объектно-ориентированных программных средств для решения комплексных задач гидродинамики на параллельных системах. -Сборник трудов Всероссийской научной конференции "Научный сервис в сети Интернет: многоядерный компьютерный мир", М.:Изд-во МГУ, 2007 г., с. 196-197.
Благодарности Работа выполнена под руководством зав. отделом ИММ РАН д.ф.-м.н., профессора Гасилова В.А., которому автор выражает искреннюю благодарность за поддержку, постоянное внимание, бесконечное терпение и бесценную помощь на всех этапах работы над диссертацией.
Автор признателен к.ф.-м.н. Карташевой Е.Л. за консультации во время работы над моделью NURBS, за помощь в разработке алгоритмов поиска прообраза NURBS кривых и поверхностей, за советы в разработке алгоритмов построения поверхностных сеток.
Автор также выражает благодарность сотрудникам ИММ РАН к.ф.-м.н. Ольховской О.Г. за помощь в разработке КО схемы теплопроводности, к.ф.-м.н. Головину М.В. за сотрудничество в разработке КЭ препроцессора, расчётной КЭ схемы диффузии, уравнений фильтрации, а также советы относительно расчётной схемы на NURBS элементах.
Автор выражает особую признательность Кузнецову Е.В. за помощь в разработке модели притока флюида к наклонно-горизонтальной скважине и решении уравнений фильтрации.
Автор благодарит сотрудников ИММ РАН к.ф.-м.н. Болдарева А.С. за помощь в реализации интерфейса произвольного граничного условия, к.ф.-м.н. Круковского А.Ю. за советы относительно метода установления, Бобкова В.Г. за огромную помощь в разработке алгоритмов STEP импорта и к.ф.-м.н. Дьяченко СВ. за моральную поддержку и многочисленные советы по разработке расчётных схем.
ПОДГОТОВКА ДАННЫХ 18
Обмен данными между системами геометрического моделирования
Высокий уровень интеграция систем CAD/CAM/CAE приводит к необходимости обмена данными технических описания продукта. Такими данными CAD системы должны обмениваться между собой или же экспортировать их в прикладные программы. Однако, все CAD системы хранят результаты проектирования в своих собственных структурах данных, формат которых зависит от конкретной системы. Они могут не соответствовать входному формату прикладной программы. Таким образом, когда две или более CAD/CAM/CAE системы объединяются в единое приложение для совместного использования данных, часто возникает проблема обмена, для решения которой необходимо иметь возможность преобразовывать данные технических требований одной системы в форму, понятную для других и наоборот. Для того чтобы облегчить преобразование и не разрабатывать программы конвертеры для всех возможных пар САПР было предложено несколько стандартных форматов обмена данными технических требований. В качестве форматов и языков межпрограммных обменов используются IGES, DXF, Express (стандарт ISO 10303-11, входит в совокупность стандартов STEP), SAT (формат ядра ACIS) и др. Такие форматы, как IGES или DXF, описывают, в основном, только геометрию объектов. Диалекты языка Express считаются более перспективными, что объясняется большей гибкостью стандартов STEP, их направленностью на различные приложения, а также на использование в современных проектных и производственных системах.
Рассмотрим кратко основные форматы обмена с точки зрения CAD систем. Форматы IGES и DXF поддерживаются всеми CAD системами и многими препроцессорами. Например, известный генератор сеток GiD 7.2 (http://gid.cimne.upc.es/) поддерживает импорт IGES и DXF, но не работает с файлами формата STEP.
Формат STEP активно используется автором при подготовке расчётных областей в рамках нескольких контекстов представления, соответствующих геометрическим моделям представления форм (раздел 1.1.1.).
Формат IGES IGES (Initial Graphics Exchange Specification) стал первым стандартным форматом обмена данными, разработанным для нужд передачи данных технических требований между различными САПР. Первая версия IGES была разработана в 1980 г. техническим комитетом под руководством Boeing, General Electric и Национального бюро стандартов США. Ранние версии IGES были неявным образом ориентированы на
CAD/CAM-системы 1970-х и начала 1980-х гг., то есть, главным образом, на обмен чертежами. В более поздних версиях спектр типов данных, подлежащих обмену, был расширен. Например, версия 2.0 поддерживает обмен данными анализа по методу конечных элементов и данными печатных плат; в версии 3.0 расширены возможности пользовательских макрокоманд, играющих важную роль при обмене стандартными библиотеками деталей; в версии 4.0 введена поддержка дерева CSG; в версии 5.0 появилась обработка данных структуры B-Rep.
IGES-файл состоит из шести разделов, которые должны идти в следующем порядке: Flag (Флаг, необязательный раздел), Start (Начало), Global (Глобальные данные), Directory Entry или DE (Запись в каталоге), Parameter Data, ли PD (Параметрические данные) и Terminate (Конец). Пять обязательных разделов идентифицируются буквами S, G, D, Р и Т.
Данные в IGES-файле могут быть представлены в двух форматах: ASCII и бинарном. Формат ASCII имеет две разновидности: фиксированную длину строки 80 символов и сжатую форму. Сжатая форма — это не что иное, как ASCII-файл, сжатый путем устранения пробелов между записями. Раздел Flag состоит из одной записи (строки) и предшествует разделу Start. В бинарном формате первый байт (8 бит) раздела Flag содержит ASCII-символ В, служащий идентификатором типа файла. В разделе Start дается описание файла в форме, воспринимаемой человеком. В нем указывается система, являющаяся источником данных, препроцессор и описываемый продукт. Записи, у которых в столбце 73 стоит символ S, принадлежат разделу Start. Раздел Global содержит информацию о препроцессоре, а также информацию, необходимую постпроцессору для интерпретации файла. В частности, в этом разделе указаны символы, используемые в качестве разделителей, формат представления целых чисел и чисел с плавающей точкой в системе-источнике, единицы измерения модели.
Раздел Directory Entry содержит список всех элементов и некоторых их атрибутов. В IGES-файле все данные технических требований представлены в виде списка предопределенных элементов: геометрических (линии, кривые плоскости, поверхности) и пояснительных (комментарии и значения размеров). Каждому элементу присваивается определенный номер типа.
Набор поддерживаемых геометрических элементов зависит от версии IGES. В частности, поддержка кривых и поверхностей NURBS есть уже в версии 2.0, а твердые тела введены в качестве геометрических элементов только в версии 4.0.
Нахождение прообраза точки В-сплайн кривой
В общем случае В-сплайн кривая описывается вектор функцией F(t) скалярного аргумента t, которая отображает отрезок в одномерном параметрическом пространстве U, t Є U в трехмерное пространство моделирования Е3, F Є Е3. Задача нахождения прообраза В-сплайн кривой на отрезке [mm; maa;] формулируется следующим образом: для заданной точки Р пространства (х, у, z) Є Е3 и точности є найти параметрический прообраз t Є U, такой что F(f) — Р \ є. Для обеспечения однозначности решения поставленной задачи мы вводим следующие ограничения на функцию F(t). Функция F(t) задаёт взаимно однозначное отображение. Тем самым, рассматриваются кривые, не содержащие самопересечений. Кроме того, F(t) - В-сплайн кривая порядка к в трехмерном евклидовом пространстве Е3, такая что кратность узловых значений не превышает к — 1 за исключением первого и последнего узла. Это обеспечивает непрерывность функции F(t). Далее, первая производная F (t) отлична от нуля. В этом случае, исключаются из рассмотрения все особые точки, например, кратные точки и точки возврата [65, 66]. При этом предполагается, что более сложные случаи обрабатываются путем предварительного подразбиения кривой (например, в точках самопересечения) на сегменты, удовлетворяющие указанным выше требованиям.
Пример, иллюстрирующий поставленную задачу поиска прообраза точки на В-сплайн кривой, показан на рис. 1.14. Изображенная кривая близка к самопересекающейся, но удовлетворяет всем вышеприведенным условиям. При этом, точка Р может не лежать строго на кривой, но находиться внутри выпуклой оболочки определяющего её многоугольника (рис. 1.15). Методы вычисления выпуклых оболочек множеств рассмотрены в книге [5].
Обозначим искомое значение параметра через t и будем искать его начальное приближение. После этого можно воспользоваться любым итерационным методом, который обеспечит определение искомого t с наперед заданной точностью є, например, методом Ньютона.
Метод поиска начального приближения для В-сплайн кривой основан на её дискретизации. В результате разбиения кривой получается конечная последовательность вершин {Pi}, соответствующая последовательности параметрических значений {U}; I = 1,..., пр разбиения кривой. Именно она используется для поиска прообраза точки Р на кривой. По построению данная последовательность включает значения параметра из узлового вектора Т = {to,...,tr}. Обозначим через Л множество индексов точек ломаной, т.е. Л = {1,...,пр}. Обозначим через d(-,-) евклидову метрику в трехмерном пространстве, в котором расположена кривая.
При этом, начальное приближение, найденное с точностью 0, можно определить как значение параметра t, соответствующее некоторой вершине Pi полученной ломаной, если выполнено условие: \F(tl)-P \ = \Pl-P \ = d(PhP ) 9. (1.12) Алгоритм поиска начального приближения включает следующие основные шаги:
1. В случае В-сплайн кривых удобно воспользоваться тем, что кривая проходит внутри выпуклой оболочки определяющего многоугольника. Поэтому вначале проверяем принадлежность точки Р указанной выпуклой оболочке с учетом 9. Если точка Р лежит вне выпуклой оболочки на расстоянии больше 9, то условие (1.12) оказывается недостижимым, в противном случае переходим к следующему шагу.
2. Вычисляем расстояния d(Pi,P ) от вершин Pi, I = 1,...,пр до точки Р\
3. Находим три ближайшие к Р точки ломаной: d{Pi, Р ) = тщ = min/eA d(Pu Р ), d(Pj, Р ) =т2 = mml=A\{i} d{Ph Р ), d{Pk, Р ) = m3 = min/=A\{ij} d(Pt, Р ).
4. Находим отрезки ломаной, инцидентные хотя бы одной из точек Pi, Pj и Pjt, и определяем минимальную из длин найденных отрезков. Обозначим её Lmjn.
5. Для всех оставшихся отрезков ломаной, длина которых больше Lmin+6, проводим оценку расстояния до точки Р . Если найденное расстояние меньше ті, то соответствующий отрезок разбивается. Если ни один отрезок не был разбит, переходим к следующему шагу. В противном случае возвращаемся к шагу 2.
6. Проверяем, являются ли Pi, Pj и Рк смежными вершинами ломаной. Если да, то переходим к следующему шагу. В противном случае, разбиваем отрезки ломаной, инцидентные хотя бы одной из указанных точек, и возвращаемся на шаг 2.
7. Проверяем условие (1.12) для точки Pj. Если оно выполняется, то Pj - найденное начальное приближение, в противном случае переходим к следующему шагу.
8. Разбиваем отрезки ломаной, инцидентные Pj, Pj и Р&, и возвращаемся на шаг 2, если не превышено максимальное число итераций алгоритма. Если число итераций больше предельного, то в качестве начального приближения выбирается Pj.
В данном алгоритме используется разбиение В-сплайн кривой и осуществляется локальное переразбиение с более мелким шагом, в случае необходимости. При этом следует задать предельное число итераций алгоритма, чтобы избежать его зацикливания в случае невозможности достижения заданной точности 9, который возникает, например, если исходная точка Р не лежит на кривой и находится на расстоянии, превышающем 9.
Классы элементов, используемые в программном комплексе
Методы Ритца и Галеркина как основа МКЭ С математической точки зрения МКЭ представляет собой обобщение метода Рэлея-Ритца-Галёркина [104], поэтому он применим к широкому классу уравнений в частных производных. Рассмотрим кратко методы Ритца и Галеркина.
Метод Ритца не производит непосредственного решения дифференциального уравнения. Вместо этого исходная задача представляется в эквивалентной вариационной формулировке, а затем ищется приближенное решение последней в виде комбинации YJC№J заданных пробных функций (pj, доставляющая минимум некоторому функционалу. При этом весовые коэффициенты Cj вычисляются из вариационного принципа, соответствующего задаче.
Пусть Н - гильбертово пространство, а Нп - его подпространство. Назовем приближенным решением по методу Ритца такую функцию и11 Є Нп, что J(un) = mm J(vn), ti» = cm. (2.1)
К построению приближенного решения можно подойти несколько иначе. Рассмотрим операторное уравнение Lu = /. (2.2)
Приближенным решением по методу Галеркина назовем элемент ип є Яп, для которого невязка Lun — f ортогональна любому элементу vn Є Нп, т.е. {Lun- f,vn) = 0, VvneHn (2.3)
Рассматривая базис пробных функций в (pi,...,(pn пространстве Я", будем искать решение ип аналогичным (2.1) образом, в виде линейной комбинации пробных функций и11 = YjfyVj- Можно показать [103], что постановки задач (2.1) и (2.3), в некотором роде, эквивалентны. Однако, уравнение (2.3) не требует положительной определенности оператора L или его симметрии и следовательно называется уравнением в слабой форме. Таким образом, уравнения Галеркина это условие не минимума, а стационарной точки [104]. Квадратичный функционал, соответствующий уравнению (2.3), выглядит следующим образом: I(v) = (Lv,v)-2(f,v). (2.4) Уравнение (2.2), являющееся точкой минимума (2.4), называется уравнением Эйлера.
Основная идея, отличающая метод конечных элементов (МКЭ) от других реализаций методов Ритца и Галеркина, заключена в том, что область, в которой требуется найти решение краевой задачи, разбивается на подобласти простой формы, называемые конечными элементами. В качестве пространства Я", в котором ищется приближенное решение, берется пространство так называемых "кусочных функций", определяемых по-своему на каждом конечном элементе и представляющих собой достаточно простые функции, например, многочлены низкой степени.
Таким образом, цель состоит в том, чтобы выбрать пробные функции достаточно удобными для вычисления и минимизации вариационного функционала и в то же время обеспечить хорошее приближение неизвестного решения.
Выбор пробных функций является существенным вопросом при реализации метода. В методе конечных элементов эти функции кусочно полиномиальны, чем объясняется его успех. Каждая функция tpj равна нулю на большей части области и отлична от нуля только в окрестности одного узла. В этой окрестности pj составлена из полиномов небольшой степени, и все вычисления становятся максимально простыми. Таким образом, конечноэлементные базисные функции должны быть 1. отличны от нуля в ограниченном множестве узлов КЭ сетки 2. и иметь минимальный носитель.
Перечисленные свойства существенны с точки зрения построения матрицы МКЭ системы, ибо в этом случае она будет разреженной, т.е. содержать много нулевых элементов.
Классы элементов, используемые в программном комплексе В рамках комплекса программ MYGEOM разработаны классы следующих элементов: 1. Линейный отрезок (2D); 2. Линейный треугольный элемент (2D); 3. Линейный прямоугольный элемент (2D); 4. Линейный треугольный элемент (3D), который используется для реализации поверхностных элементов трехмерных тел, дискретизованных призматическими и тетраэдральными элементами;
5. Линейный прямоугольный элемент (3D), который используется для реализации поверхностных элементов трехмерных тел, дискретизованных кубическими и тетраэдральными элементами; 6. Линейный тетраэдральный элемент; 7. Линейная треугольная призма; 8. Линейный кубический элемент или прямоугольная призма; 9. NURBS кривая (2D);
10. NURBS поверхность (2D). NURBS элементы, реализованные в программном комплексе, созданы на основе NURBS кривых и поверхностей (см. главу 1). Вопросы подготовки расчётной NURBS геометрии описаны в разделе 2.1.5. Базисные функции лагранжевых элементов и основы МКЭ подробно рассмотрены в литературе, например, [105-107]. Вопросы, связанные с использованием граничных элементов, хорошо изложены в [108,109].
Нахождение прообраза NURBS поверхности
Рассмотрим несколько примеров работы алгоритма поиска прообраза В-сплайн кривой. Отметим, что с целью адекватного тестирования и иллюстрации были выбраны неравномерные начальные разбиения соответствующих кривых.
На рис. 3.1, 3.2 изображен достаточно простой случай, когда точка Р точно лежит на кривой. На первой итерации алгоритма на исходной ломаной найдены три последовательные точки Pj, Pj и Рк, соответствующие расстояниям mi, тг и т3 от Р . Проведено дополнительное разбиение инцидентных отрезков, поскольку основной критерий (1.12) оказался не выполненным. Алгоритм заканчивает работу на второй итерации (рис. 3.2).
На рис. 3.3, 3.4 изображен более сложный случай. Точка Р не лежит на кривой, но находится в пределах в от неё внутри выпуклой оболочки вершин определяющего многоугольника. Начальное приближение найдено за несколько итераций алгоритма. При этом, были дополнительно разбиты только отрезки, инцидентные точкам PitPj и Рк, найденным на первой итерации. Дополнительного разбиения других отрезков ломаной не потребовалось.
Самый сложный случай изображен на рис.3.5, 3.6. На некотором этапе работы алгоритма были найдены три последовательные точки ломаной, соответствующие mi, тг и тз- Далее, производилась оценка расстояний до всех отрезков текущей ломаной, длина которых больше Lmin-\-6, в результате возникла необходимость разбиения отрезка ДРг Как видно из приведенных примеров, скорость и точность нахождения начального приближения зависят от исходного разбиения кривой. Для обеспечения адекватного решения необходимо, чтобы точность аппроксимации кривой ломаной была не ниже заданной точности искомого начального приближения. При этом, если выбрать в достаточно маленьким, т.е. близким к є, то данный алгоритм дает решение задачи поиска прообраза и нет необходимости дополнительно использовать другие итерационные методы.
Нахождение прообраза NURBS поверхности Рассмотрим работу алгоритма для рациональной В-сплайн поверхности, изображенной на рис. 3.7. Данная поверхность определяется следующим набором параметров: пи = 3 - верхний индекс числа контрольных точек в направлении и, nv = 3- верхний индекс числа контрольных точек в направлении v, к = 4 - порядок кривой в направлении и, 1 = 3- порядок кривой в направлении v, (О, 0, 0, 0, 1, 1, 1, 1) - узловой вектор в направлении и, (О, 0, 0, 1, 2, 2, 2) - узловой вектор в направлении v,
Выберем точки, которые лежат на некотором расстоянии от поверхности, что наиболее типично для практических задач. Точка Pi(13.5,0.7125,-15.9975) лежит вблизи границы поверхности, точка Рг(—13,2.35484,5.56452) находится над поверхностью, точки Р3(-15.1,4.55413,3.13394) и Р4(-15.1,4.66413,3.13394) находятся вблизи изгиба поверхности. Во всех примерах начальный шаг разбиения по и равен 0.1, а по v - 0.2, а точность поиска начального приближения в = 0.5.
В результате работы алгоритма для точки Pi получено начальное приближение:(u,v) = (0.95,1.9), что соответствует минимальному расстоянию 1.008. Для точки Рг найдено начальное приближение (и, v) = (0.05,0.48), что соответствует минимальному расстоянию 2.41925. Таким образом, в первых двух примерах заданная точность 9 не достигнута. Полученный результат обусловлен тем, что указанные точки находятся слишком далеко от поверхности.
Для точек Рз и Рь найдены начальные приближения (0.6,1.84) и (0.05, 0.81) соответственно. Минимальные расстояния до поверхности составляют 0.347 в первом случае и 0.278 - во втором, что меньше заданной точности в. Результат работы алгоритма для двух последних примеров приведён на рис.3.9.
