Содержание к диссертации
Введение
1 Компьютерная алгебра и разработка программных комплексов численного моделирования 14
1.1 Применение систем компьютерной алгебры в задачах МКЭ 14
1.2 Разработка интерфейс-преобразователя 16
1.3 Метод взвешенных невязок и альтернативные интегральные формы 18
1.4 Примеры постановки и решения сложно-сопряженных задач 26
2 Автоматическое порождение текстов программ 45
2.1 Основные определяющие соотношения, модели сред 45
2.2 Функциональные зависимости теории пластичности 50
2.3 Схема построения тангенциальной матрицы жесткости 57
2.4 Стратегия символьных преобразований 60
2.5 Замечания по процедурам символьной генерации программ. 70
3 Объектно-ориентированное структурирование МКЭ 76
3.1 Требования к программным комплексам МКЭ 76
3.1.1 Вычислительные требования 76
3.1.2 Требования к структуре программы 80
3.1.3 Требования к виду шаблона системы МКЭ з
3.2 Представление структуры метода конечных элементов 83
3.3 Объектно-ориентированное программирование МКЭ 90
3.4 Классы в МКЭ 97
3.4.1 Структура классов 98
3.5 Требования, предъявляемые к OOFEM 103
3.6 Классы метода конечных элементов.
3.6.1 Методы 109
3.6.2 Класс Node и его основные методы. 111
3.6.3 Класс Element и его основные методы 118
3.6.4 Класс Material и его основные методы 125
3.6.5 Класс (свойств) Property и его основные методы 126
3.7 Схема хранения модели 126
3.7.1 Класс List (связанные списки) 127
3.7.2 Применение связанных списков
3.8 Алгебраические классы 129
3.9 Пример применения шаблона
3.9.1 Ввод модели 134
3.9.2 Генерация модели 136
3.9.3 Формирование глобальной системы уравнений 137
3.9.4 Решение глобальной системы уравнений 138
3.9.5 Обработка результатов 139
3.10 Линейные и нелинейные модели. Динамические задачи. 142
Основные результаты и рекомендации. 145
Литература
- Разработка интерфейс-преобразователя
- Примеры постановки и решения сложно-сопряженных задач
- Схема построения тангенциальной матрицы жесткости
- Объектно-ориентированное программирование МКЭ
Разработка интерфейс-преобразователя
Наличие ошибок, возникающих в процессе программирования, общеизвестно. В среднем каждый программист совершает не менее двух существенных ошибок при написании одной страницы кода программы в соответствии с заданной спецификацией, или, используя другой критерий эффективности, средний программист пишет не боле 10-15 отлаженных операторов в день, что указывает на крайнюю неэффективность работы. Научное сообщество пыталось исправить это путем введения различных парадигм программирования, таких как структурное или объектно-ориентированное. Однако создание программных продуктов в наукоемких областях сильно осложняется необходимостью не только отображать сложные логические связи спецификаций решаемых задач, но и преобразовывать математические модели, записанные на своих собственных языках, наиболее удобных для их представления в языки программирования, возникшие на базе других предметных областей. К языкам представления моделей для задач механики сплошных сред можно отнести тензорное исчисление и дифференциальные уравнения в частных производных. Очевидным желанием людей, связанных с разработкой программных продуктов является автоматизация процесса программирования, особенно в тех обла 15
стях, где создание прототипов занимает существенную часть работы. Это особенно важно при программировании вычислительных моделей, поскольку вычислительная эффективность решения конкретной задачи, например, при численном решении уравнений в частных производных зависит от применяемой схемы дискретизации, используемой библиотеки методов вычислительной линейной алгебры, формы представления данных, т.е. схем хранения разреженных матриц и т.д. Не последним является и вопрос о возможности эффективного распараллеливания построенной программы. В этом разделе в первом приближении будут рассмотрены некоторые вопросы использования систем компьютерной алгебры (СКА) для порождения программ МКЭ. Здесь не будут затрагиваться вопросы построение библиотек различных конечных элементов и процедура дискретизации области решения, хотя применение СКА в этих областях достаточно прозрачно. Можно указать на три области потенциального применения СКА в МКЭ:
В качестве интеллектуального интерфейса, позволяющего преобразовывать входной язык математической модели, например, дифференциальных уравнений в выходную последовательность утверждений, представляющих собой программные модули вычислительной модели, представленные на одном из языков программирования высокого уровня. В данной работе к качестве такого языка выбран Фортран 95. Хотя можно использовать и любой другой язык, содержащий достаточное, для решения задачи, подмножество средств объектно-ориентированного программирования.
Совмещение операторов вычислительно языка типа Фортран с one 16 раторами символьных преобразований в один текстовый модуль для автоматической генерации недостающих частей программ в вычислительном комплексе, что будет проиллюстрировано на примере генерации модуля соответствующего тангенциальной матрице жесткости для части вычислительной модели, описывающей упругопласти-ческие деформации скелета горной породы.
Автоматизированная генерация матричного представления МКЭ. Последняя проблема достаточно давно и хорошо описана в научной литературе /32, 44/ и др., поэтому здесь она рассматриваться не будет.
Определим первую задачу автоматизации разработки программных комплексов численного моделирования следующим образом: необходимо построить интерфейс для представления спецификации исходной задачи в форме близкой к языку предметной области. В нашем примере это некоторая система дифференциальных уравнений в частных производных. Разработанный пакет представляет собой набор процедур, написанных на языке системы Maple и модулей Фортран 95, которые используются для автоматической генерации программ численного решения систем уравнений в частных производных методом конечных элементов. Пакет предназначен для решения эволюционных систем дифференциальных уравнений методом расщепления по процессам. В случае, если уравнение записывается для векторных или тензорных величин система решается покомпонентно, что также можно рассматривать в качестве схемы расщепления. Хотя в данной работе вопросы сходимости и устойчивости не рассматриваются, но поскольку любую схему расщепления можно интерпретировать, как некоторую форму блочно-итерационного метода решения систем линейных или нелинейных алгебраических уравнений, то ясно, что анализ качества алгоритма можно осуществить, опираясь на указанное представление. Пространственную дискретизацию области решения задачи можно провести с помощью любой программы построения неструктурированных сеток. В этом случае необходимо написать соответствующий ей модуль ввода геометрических координат узлов элементов и их инцидентности. В пакете структура сетки узлов конечных элементов представлена в форме связанных списков, что позволяет проводить динамическое реструктурирование конечно-элементного разбиения в процессе счета, например, можно добавить элементы в области больших градиентов решения. В настоящее время в пилотной версии пакета используется некоторая стандартная библиотека одномерных, двумерных и трехмерных элементов: одномерный линейный элемент, линейный треугольный элемент и линейный тетраэдр для решения трехмерных задач. Вся работа по генерации кода программы, компиляции, решению и визуализации полученных результатов осуществляется автоматически начиная с вариационной формы представления МКЭ, которая задается пользователем и, соответственно, представляет существенную часть спецификации решаемой задачи. МКЭ, в данной работе, рассматривается, в качестве специальной формы метода взвешенных невязок, что позволяет без существенных изменений обобщить данный подход на другие схемы решения типа метода граничных элементов, различного рода кол локаций и др. Таким образом, спецификация задачи состоит в определении системы решаемых уравнений, вводе граничных условий и указании программы решения, полученной в результате дискретизации, системы алгебраических уравнений. В текущей версии пакета используются условия Дирихле, Неймана и Робена (смешанные). Можно расширить пакет для применения граничных условий, определенных в /4/, которые необходимы для корректного численного моделирования волновых процессов.
В общем случае система преобразований должна быть, построена так, чтобы предоставить пользователю максимальные возможности влияния на построение вычислительной схемы. Поэтому ввод решаемой системы уравнений состоит в представлении ее в эквивалентной вариационной форме. Здесь эквивалентность понимается с точностью до соответствия слабого или обобщенного решения классическому. Этот подход можно определить следующим образом.
Примеры постановки и решения сложно-сопряженных задач
Далее для твердой фазы необходимо выписать определяющее соотношение, если модель предназначена для изучения пластических деформаций скелета горной породы необходимо использовать инкрементальную теорию, связи прироста напряжения с приростом деформации, если модель предназначена для описания упругих (даже нелинейных) деформаций матрицы пористой среды достаточно использовать деформационную теорию связи напряжения с деформацией, т.е.
В первом случае матричное представление тензора D называется танген-цальной матрицей жесткости. Аналогичным образом вводится связь между вектором смещений и тензором деформаций, например, для случая малых деформаций имеем:
Если рассматривать обе задачи, первую в предположении пренебре-жимой малости VT, а вторую для упругого взаимодействия, полагая, что горизонтальные смещения и скорости изменения деформации пренебрежимо малы, то в результате получится система эллиптических и параболических уравнений.
В качестве демонстрационного примера, иллюстрирующего работу интерфейса-преобразователя, рассмотрим решение следующей системы уравнений в частных производных:
Для определенности будем предполагать, что S\ - граница на которой заданы условия Неймана, a, S2 граница с условиями Дирихле. Рассмотрим пример решения системы дифференциальных уравнений в частных производных эллиптического и параболического типа, соответственно: где индекс і означает соответствующий временной слой, a At - шаг по времени. Отметим, что для аппроксимации производных по времени можно использовать любые схемы и шаблоны, принятые в численных методах решения обыкновенных дифференциальных уравнений.
Интеграл по границе, с учетом того, что S = S\ \J 5г и S\ - граница на которой заданы условия Неймана, а 5г - граница с условиями Дирихле, можно представить в виде s ф %dS — s ф dS + 3ф 5. Полагая, что весовая функция ф = 0 на S2, получим $s ф ;dS = s ф %;dS. Именно эта форма интеграла по границе и будет использоваться в дальнейшем.
Необходимо отметить, что для простоты изложения, схема построения вычислительной модели была ограничена только методом Галеркина с интегрированием по частям. В пакете реализованы аналогичные вычислительные схемы, основанные на различных формах метода коллокаций, метода наименьших квадратов и метода граничных интегральных уравнений в форме граничных элементов. Вид схемы собственно зависит от выбора систем весовых и пробных функций.
Кроме того в пакете реализована, описанная выше, процедура метода множителей Лагранжа, что позволяет автоматически строить различные представления системы дифференциальных уравнений для одной и той же физической задачи, путем построения соответствующих функционалов исходя из уравнений сохранения, равновесия и определяющих соотношений. Например, можно различным образом объединять уравнения сохранения массы с уравнением Дарси или другим определяющим соотношением, связывающим скорость фильтрации с давлением и получать различные представления уравнений фильтрации.
В результате выполнения указанных действий для построения вычислительной модели необходимо ввести следующую часть рабочей страницы в систему Maple:
В данном примере независимые переменные и и v определены явно, как величины зависящие от пространственных координат т.е. и(х,у), v(x,y). Значение v на предыдущем временном шаге обозначается через vO(x, у), ф(х,у) - весовая функция.
Ключевые слова (операторы) E_Int и B_Int используются для представления интегралов по области решения (элементу) и по ее границе, что указывает на необходимость их вычисления при обращении к основной процедуре обработки входных строк. Все величины, относящиеся к описанию уравнения должны находиться в области действия операторов E_Int и B_Int.
Запись исходных уравнений можно значительно упростить, используя инфиксный оператор Szts, реализующий соглашение о суммировании при вычислении произведения декартовых тензоров, оператор nab, соответствующий V, и известные тензоры delta[i,j] = S j и epsilon[i, j, к) = ei,j,k- В результате пример примет следующий вид:
Схема построения тангенциальной матрицы жесткости
Разработка новых определяющих соотношений (моделей) используемых, в частности, для разработки новых конечно-элементных систем моделирования представляет очень важную область исследования в технических дисциплинах. Это подтверждается интенсивными исследованиями, ведущимися в этом направлении в таких технических областях, как жаростойкие композиты, железобетонные конструкции, геотехнические среды /45/, к которым относятся коллекторы нефти и газа, и другие. Исследования в области определяющих соотношений включают в себя разработку математических моделей для предсказания нелинейной реакции сред на внешние воздействия. Обычно это осуществляется путем введения различного рода тензоров (матриц), связывающих тензор деформации с тензором напряжения в форме обобщенного закона Гука. Матричная форма представления тензорных величин более удобна для реализации вычислений методом конечных элементов. Например, для инкрементальных моделей используют тангенциальные матрицы жесткости или податливости среды.
Далее на основе математической модели разрабатывается спецификация вычислительной модели (алгоритм), которая служит основой разработки программного кода. Финальной частью проекта является проверка кода, т.е. соответствия программы ее спецификации. Последняя задача, в общем случае, исключительно трудна для разрешения, если не использовать для написания программ языков, позволяющих отразить внутреннюю логику спецификаций, например, Пролог, Схема, РОР-11 и др. или применить систему автоматической генерации текста программы.
Очевидно, что полный цикл построения программного обеспечения требует проведения большого количества рутинных алгебраических операций и программирования. Поэтому время, затрачиваемое на получение результата в соответствии с указанной последовательностью действий большое, порядка многих месяцев. Из этого следует, что для исследователя трудно внести кардинальные изменения в теорию, задающую определяющие соотношения, так как для этого требуются существенные усилия. Кроме того, результат работы подвержен естественным человеческим ошибками, которые часто трудно обнаружить. В силу указанных причин, ясно, что только компьютерные аналитические преобразования могут обеспечивать существенный стимул к развитию и реализации различных теорий определяющих соотношений и их конечно-элементных реализаций. Для определенности рассмотрим в качестве примера построение тангенциальных матриц жесткости в рамках инкрементальной теории упругопластических деформаций, хотя очевидно, что данный подход можно использовать и для инкорпорирования других определяющих соотношений в вычислительные модели.
Язык СКА MAPLE дает возможность эффективным способом получать характеризующие среду матрицы для новой определяющей модели. Существует несколько очевидных преимуществ использования СКА: существенное уменьшение рутинной работы при генерации тангенциальной матрицы жесткости, увеличение надежности полученных результатов анализа и небольшое время создания новой модели определяющих соотношений. Однако, прямое, лобовое использование любой СКА не возможно. Главное препятствие - экспоненциальный рост алгебраических выражений при выполнении промежуточных действий, что требует существенных объемов памяти и времени счета. Кроме того, возникают проблемы, связанные с автоматическим переводом полученных соотношений в программу на одном из процедурных языков программирования. Например, непосредственное использование системы MAPLE для построения тангенциальной матрицы жесткости трехмерной модели упругопластических деформаций Мора-Кулона с оптимизацией кода на языке Фортран-77 приводит к программе, состоящей примерно из 5000 операторов, что абсолютно недопустимо.
Символьных преобразования для построения программ МКЭ применялись и ранее. Однако большинство работ сконцентрировано на автоматическом порождении матриц инерции (масс) и матриц жесткости /32, 40, 42, 43, 44/. Последние рассматриваются в смысле стандартной терминологии МКЭ.
В данной части работы будет продемонстрировано: применение СКА MAPLE для порождения определяющих матриц среды, стратегия, позволяющая упрощать промежуточные выражения, и генерация модуля программы на языке Фортран-95, ответственного за соответствующие вычисления. Разработанная процедура использовалась для порождения тангенциальных матриц жесткости для следующих моделей упругопластических деформаций: пластичность металла по Мизесу, модель пластичности грунта Друкера-Прагера, модели пластичности бетона и других твердых сред Они интенсивно применялись в для исследований упругопластического поведения различных сред и в геотехнических приложениях. Необходимо отметить, что в процессе тестирования были обнаружены ошибки в опубликованных формулах, которые далее перекочевали в программные продукты, например, в широко используемой работе /33/. Ниже приведены различные формы критерия текучести в виде некоторой функции тензора напряжений, заданной в некотором пространстве инвариантов для моделей на которые применялись для тестирования разработанной процедуры.
Объектно-ориентированное программирование МКЭ
Концепция класса ООП дает возможность использовать символическое программирование, применяя перегруженные операторы. Перегруженный оператор - это переопределение обычного оператора для адаптации его к текущему контексту. Некоторые типы переменных естественным образом подходят для построения перегруженных операторов. Это особенно относится к математическим типам, таким как рациональные числа, комплексные числа, векторы, матрицы, тензоры и т.д. В этих случаях имеет смысл использовать общеупотребительную инфиксную запись бинарных операторов, которая уже перегружена даже для простых типов переменных, какими являются целые числа, вещественные числа и т.д. Перегрузка оператора делает возможным установить прямую связь между математическими понятиями и применением их в естественно-научном программировании. Классы типа Матрица и Вектор могут использоваться как инструментальные средства во всех парадигмах программирования. По существу они не являются иерархическими классами в общепринятом объектно-ориентированном смысле, например, класс IntegerVector не может быть получен из RealVector. Такие классы названы Б.Страуструпом абстрактными типами данных (АТД) /52/.
В некоторых работах рассматривается возможность использования Матричных и Векторных классов в программах на языке C++. Большинство авторов применяет их как средство программирования в рамках полностью объектно-ориентированной структуры шаблона /35, 50, 54/. Другие авторы изучали влияние, которое эти классы могут оказать на разработку процедурного, научно-технического программного обеспечения /48, 49/. Везде в качестве вывода отмечалось, что использование алгебраических инструментальных средств увеличивает эффективность разработки и делает код менее подверженным ошибкам. Также, указано на улучшение удобочитаемости программ благодаря более символическому (алгебраическому) стилю программирования.
Предшествующую часть работы можно использовать в качестве базиса для построения объектно-ориентированного шаблона МКЭ, который мы назовем OOFEM. Сначала необходимо определить структуру простой программы, например, в рамках традиционной задачи линейной теории упругости и требования которым должна удовлетворять, структура управления. Линейная теория упругости рассматривается в качестве наиболее простой, но содержащей все необходимые черты предметной области в которой моделирование на основе МКЭ доказало свою состоятельность и эффективность.
Введем следующие основные классы МКЭ: Node (Узел), Element (Элемент), Material (Материал) и Property (Свойства). Ниже будет приведено описание их функционирования в шаблоне ООП, их основные атрибуты и методы. Будут определены классы List и алгебраические классы (Matrix) Матрица и (Vector) Вектор. В заключение, в качестве примера, будет представлен псевдокод для задачи линейной теории упругости и указано на те части модели, которые необходимо изменить при введении новых элементов, материалов или решении нелинейных задач, например, при моделировании нелинейных сред . Как было указано ранее, в МКЭ можно выделить четыре основных класса в ООП шаблоне OOFEM: Node, Element, Material и Property.
Класс Node играет двойную роль с одной стороны это точка сопряжения элементов (геометрический узел), с другой это интерполяционный узел, отвечающий за работу с неизвестными, решаемой задачи, называемыми степенями свободы - dof и нагрузками - load. В дальнейшем, для простоты, будем полагать совпадение этих типов узлов. Для того, чтобы в глобальной системе число уравнений равнялось числу неизвестных необходимо в каждом узле задать один элемент из дуальной пары (dof, load). Эта дуальность влияет на решение глобальной системы уравнений, потому что задание dof вводит неизвестную начальную нагрузку, поэтому глобальная переменная dof и глобальный вектор нагрузки одновременно могут содержать неизвестные значения. Предлагается указанную дуальность учитывать на уровне методов решения так, чтобы метод Matrix_solve одновременно вычислял неизвестные составляющие степеней свободы и реакций. Node должен передавать заданные dof и load в глобальную систему уравнений и доопределять неизвестные значения, полученные в результате ее решения. Node также определяет точку, принадлежащую области решения задачи, задаваемую вектором абсолютных координат (глобальные геометрические координаты узлов). Это различие между узлами сопряжения и геометрическими узлами можно использовать для введения подчиненных, промежуточных узлов элементов.
