Если у вас возникли сложности с курсовой, контрольной, дипломной, рефератом, отчетом по практике, научно-исследовательской и любой другой работой - мы готовы помочь.
Предоплата всего
Подписываем
МАТЕМАТИЧЕСКОЕ ОБЕСПЕЧЕНИЕ САПР
Математическое и программное обеспечения САПр призваны решить две задачи.
В соответствии с этим математическое обеспечение САПр рассмотрим применительно к модулям CAD и CAE.
МАТЕМАТИЧЕСКОЕ ОБЕСПЕЧЕНИЕ МОДУЛЯ CAD
Эффективность и удобство создания модели определяются возможностями математического и программного обеспечения, поддерживающего модуль CAD. Одна из главных функций этого программного обеспечения – создание графического изображения на экране дисплея.
Принципы формирования изображения
Геометрические задачи
При формировании изображения в первую очередь решаются две геометрические задачи:
На рис. 1 показаны следующие системы координат:
С каждым элементом модели (стержнем, поверхностью, опорой и т. д.) связана своя локальная система координат. Локальные системы координат вводятся для удобства программирования процедур, осуществляющих прочностные расчеты, и сейчас не рассматриваются.
Работать с моделью удобно в системе координат XoYoZo. Получить проекцию модели на экране дисплея удобнее в том случае, если координаты характерных точек модели вычислены в мировой системе координат.
Рассмотрим алгоритм вычисления координат точки в мировой системе координат на примере плоской модели (рис. 2).
Расположение системы координат модели (XoYo) в мировой системе координат (XwYw) описывается тремя параметрами:
Точка в системе координат модели имеет координаты Xo, Yo.
Требуется вычислить координаты Xw, Yw.
Как следует из рис. 2,
Xw = Xc + Xo cos - Yo sin
Yw = Yc + Xo sin + Yo cos
или, в матричной форме
.
Матрица преобразования формируется по трем параметрам Xc, Yc, и применяется ко всем характерным точкам модели.
Известно, что поворот вокруг любой произвольной оси раскладывается на повороты вокруг осей X, Y, Z. Таким образом, расчет координат при повороте пространственной модели вокруг произвольной оси можно представить совокупностью рассмотренных выше элементарных операций.
После вычисления координат точек модели в системе XwYwZw можно формировать изображение на экране дисплея.
Рассмотрим алгоритм вычисления координат U, V на экране дисплея какой-либо точки 1, принадлежащей модели (рис. 3).
В мировой системе координат определены следующие данные:
Требуется определить координаты X1, Y1, Z1 в мировой системе координат и координаты точки U1, V1 в системе UV.
Из подобия треугольников следует: . Из этого уравнения вычисляется X1. Соответствующая координата U1 в системе координат дисплея находится по формуле U1 = X1 – Xd.
Аналогично вычисляется координата V1.
Задача о видимости точек
Далее при необходимости решается задача о видимости точек.
Существует множество алгоритмов получения изображения модели со скрытыми невидимыми элементами.
Рассмотрим «Алгоритм художника».
Элементы модели (поверхности, линии, точки) сортируются по удаленности от точки зрения. Рисование начинается с самого дальнего элемента, который имеет точку с максимальной координатой Yw (рис. 1). Если ближние элементы проецируются в ту же область экрана, что и дальние, то ближние элементы автоматически закроют дальние (частично или полностью).
При выводе сложных поверхностей используется прием, называемый триангуляцией или фасетированием. Поверхность представляется набором небольших треугольников. Каждый треугольник рассматривается как отдельный элемент, к которому можно применить алгоритм художника.
Задача тонирования
Тонирование (rendering) – это воспроизведение на экране дисплея эффектов, создаваемых светом.
Существует две основные технологии тонирования – затушевывание и трассировка лучей.
Рассмотрим технологию затушевывания.
Затушевывание аналогично удалению невидимых элементов (может применяться тот же алгоритм художника), но цвет каждой точки модели зависит от интенсивности отраженного света.
На рис. 4 показано освещение поверхности и отражение поверхностью падающего света.
Прямой свет – это свет, исходящий непосредственно от источника. Предполагается, что источник света расположен достаточно далеко (например, солнце), и падающие лучи можно считать параллельными.
Окружающее освещение – это свет, отраженный от множества других поверхностей. Предполагается, что он приходит с бесчисленного множества направлений.
В технологии затушевывания предполагается также, что наблюдатель находится далеко от объекта, и для любой точки поверхности направление взгляда одинаковое.
Отражение окружающего освещения R1 считается равномерным во всех направлениях и вычисляется по формуле
,
где I1 – интенсивность окружающего освещения, Kd – коэффициент отражения поверхности.
Наблюдатель будет воспринимать одну и ту же интенсивность отраженного окружающего освещения независимо от своего положения.
Отражение прямого света рассматривается как комбинация двух видов отражения: диффузного и зеркального.
Диффузное отражение состоит в том, что поверхность поглощает свет, а затем излучает его равномерно во всех направлениях. Интенсивность диффузного отражения R2 записывается в соответствии с «законом косинусов Ламберта»:
,
где I2 – интенсивность прямого света. Угол называется углом падения.
Зеркальное отражение – это прямое отражение света поверхностью. Существует много моделей прямого отражения. Наиболее распространена «модель Фонга», в соответствии с которой интенсивность зеркального отражения R3 записывается в виде:
,
где Ks – коэффициент зеркального отражения, f – угол между отраженными лучами и направлением взгляда.
Показатель степени n > 1 зависит от материала поверхности. Например, для металлов n > 100, а для асфальта n 1.
Интенсивность отражения в данной точке поверхности равна сумме интенсивностей: R = R1 + R2 + R3.
Для учета оттенка отраженного света сумма интенсивностей R = Ri записывается отдельно для каждого базового цвета (красного, зеленого, синего).
Каждый материал по-разному отражает базовые цвета. В свою очередь, источники света могут иметь различную интенсивность базовых цветов. Поэтому величины I1, I2, Kd зависят от базового цвета.
Зеркальное отражение практически не влияет на цвет, вследствие чего коэффициент зеркального отражения Ks и показатель степени n не зависят от базового цвета. Цветовой оттенок зеркального отражения определяется комбинацией интенсивностей I2 базовых цветов источника света.
Технология затушевывания применяется в относительно простых условиях, когда один объект освещается далеко расположенными источниками света.
В сложных ситуациях (наличия многих объектов, особенно если некоторые из них пропускают и преломляют лучи, а также близко расположенных источниках света, как, например, при тонировании интерьеров) используется технология трассировки лучей.
Рассмотрим технологию трассировки лучей.
Источники света испускают лучи во всех направлениях. Лучи прослеживаются до тех пор, пока они не пересекают плоскость дисплея. В процессе движения луча он попадает на объекты модели и изменяет свою интенсивность и направление по законам рассеяния, отражения и преломления. Но так как количество лучей равно бесконечности, такой прямой подход к трассировке затруднителен.
При практической реализации рассматриваются траектории лучей в обратном направлении от точки зрения через каждый пиксел экрана. Движение по лучу заканчивается в один из следующих моментов:
Оттенки и интенсивности отражаемого света на рассеивающих поверхностях должны быть вычислены заранее по методу затушевывания.
Очевидно, что вышеприведенные алгоритмы просты для реализации, однако могут требовать большого времени для выполнения, особенно в программах, реализующих твердотельное моделирование.
Системы геометрического моделирования
Системы геометрического моделирования делятся на следующие классы:
Термины в скобках обозначают ключевые слова, которые встречаются в названиях программных комплексов, в которых реализована та или иная система.
Системы каркасного моделирования
Каркасные модели полностью описываются в терминах точек (узлов) и линий (элементов). С каждым элементом должны быть связаны численные характеристики материала и характеристики поперечного сечения.
Существуют программные комплексы, поддерживающие моделирование криволинейных элементов переменного сечения.
Если программный комплекс позволяет описывать элементы только постоянного сечения, то такие системы моделирования в некоторой литературе относят к двухсполовиномерной (2,5-мерной) геометрии.
Каркасным моделям присуща неоднозначность в том смысле, что при отображении модели на экране дисплея возможна двоякая интерпретация изображения. Интерпретация зависит от того, какие узлы и элементы воспринимаются на переднем плане, а какие – на заднем.
Системы поверхностного моделирования
Поверхностные модели описываются в терминах точек (узлов) и поверхностей (элементов). Для элементов поверхностей необходимо задавать численные характеристики материала и толщину.
Существует три стандартных метода создания поверхностей:
Методы создания поверхностей аналогичны методам создания объемных тел в системах твердотельного моделирования, за исключением того, что создается поверхность, не имеющая толщины, а не заполненный объем.
Системы твердотельного моделирования
СТМ являются полномасштабными системами и используются, как правило, на предприятиях машиностроения.
В СТМ реализованы следующие возможности моделирования.
На рис. 5 показано выполнение булевых операций над двумя примитивами – параллелепипедом и конусом. В результате получаются новые тела нестандартных форм.
Если шаблон – фигура не замкнутая, то получится не объемное тело, а поверхность. Такой вариант реализован в системах поверхностного моделирования.
Если к натянутой поверхности не добавить две торцевые поверхности, то получится не объемное тело, а только натянутая поверхность. Такой вариант реализован в системах поверхностного моделирования.
Моделирование границ имеет другое название – «функции тонкого редактирования».
Например, изображенное на рис. 5 конусное отверстие может быть добавлено в деталь путем выбора соответствующей опции меню программы и описания параметров отверстия. Создание конуса и выполнение Булева вычитания более трудоемко.
Например, форму и положение прямоугольника на плоскости можно определить пятью параметрами, показанными на рис. 12.
Вследствие этого параметрическое моделирование получило название «моделирование по размерам».
Декомпозиционное моделирование можно сравнить с разбивкой объекта на конечные элементы, выполняемой на стадии создания формы.
К типичным декомпозиционным моделям относятся воксельное представление, представление октантного дерева и ячеечное представление.
Воксельное представление
Воксел – это трехмерный аналог пиксела. Термин получен комбинацией слогов: volume (объем) + pixel (пиксел) = voxel.
Воксельное представление – это трехмерный аналог растрового представления плоской фигуры. Объем, занимаемый моделью, делится на маленькие кубы плоскостями, перпендикулярными осям x, y, z. Расстояние между плоскостями определяется желаемой точностью воксельного представления. Модель представляется в виде трехмерного массива, количество элементов которого совпадает с количеством кубиков. Каждый элемент массива принимает значение 0 или 1 в зависимости от того, охватывает кубик область, занятую материалом модели или нет.
На рис. 15 показано довольно грубое воксельное представление конуса диаметром основания 10 единиц, высотой 8 единиц с помощью массива размерностью 10 10 8.
Воксельное представление имеет ряд преимуществ по сравнению с ранее рассмотренными методами моделирования:
Представление октантного дерева
Октантное представление аналогично воксельному представлению в том смысле, что модель представляется совокупностью кубов. Однако в октантном представлении размеры кубов переменны из-за другой схемы деления пространства.
Сначала создается куб (корневой октант), в который модель помещается целиком. Корневой октант делится на восемь октантов (2 2 2), после чего анализируется их положение относительно модели. Если октант находится полностью внутри модели или полностью вне модели, то их деление на этом заканчивается. Над другими октантами процедура деления повторяется. Деление завершается при достижении заданного минимального размера октанта.
Для хранения в памяти компьютера таких моделей используется иерархическая структура (рис. 16), что и дало название «октантное дерево».
Ячеечное представление
Ячеечное представление – это представление модели в виде комбинации более простых базовых элементов. На форму элементов не накладывается жестких ограничений.
Немногообразные системы моделирования
В математике существует термин «условия многообразности». Этим условиям отвечают модели, созданные в каркасных, поверхностных и твердотельных моделях. Нарушениями условий многообразности являются, например:
и т.д.
Немногообразные системы моделирования позволяют не соблюдать условия многообразности, а также использовать элементы каркасных, поверхностных и твердотельных моделей в одной среде проектирования.
Системы моделирования устройств
Каркасные (в меньшей степени), поверхностные и твердотельные (в большей степени) системы используются главным образом для разработки отдельных деталей.
Системы моделирования агрегатов предоставляют средства для группировки деталей в устройства и разработки соединений деталей.
Чаще всего агрегат представляется в виде дерева, которое показывает вхождение деталей в модули устройства (аналогично «Проводнику» (Explorer) операционной оболочки Windows).
Системы моделирования агрегатов являются интеллектуальными в том смысле, что изменение параметров одной детали автоматически приводит к необходимым изменениям во всем устройстве.
После сборки устройства система моделирования предоставляет возможность выполнить, например, анализ кинематики или расчет методом конечных элементов.
Информация в системах твердотельного моделирования
Информация о модели хранится в структурах трех типов.
Пример дерева C-Rep показан на рис. 13. Операция P or K соответствует рис. 5. После этого из полученного тела вычитается цилиндр C для образования отверстия.
Каждый узел дерева содержит указатели на другие узлы и на области памяти, где записана информация о размерах примитива и его расположении в пространстве, либо информация об операции над примитивами.
Для примера на рис. 14 показан тетраэдр, для описания границ которого можно пользоваться информацией, приведенной в табл. 1.
Таблица 1 – Информация B-Rep
|
Вершины |
Ребра |
Грани |
||||||||||||
|
V1 |
x1 |
y1 |
z1 |
R1 |
V1 |
V2 |
G1 |
R1 |
R3 |
R2 |
||||
|
V2 |
x2 |
y2 |
z2 |
R2 |
V2 |
V3 |
G2 |
R1 |
R5 |
R4 |
||||
|
V3 |
x3 |
y3 |
z3 |
R3 |
V1 |
V3 |
G3 |
R2 |
R6 |
R5 |
||||
|
V4 |
x4 |
y4 |
z4 |
R4 |
V1 |
V4 |
G4 |
R3 |
R4 |
R6 |
||||
|
R5 |
V2 |
V4 |
||||||||||||
|
R6 |
V3 |
V4 |
||||||||||||
Последовательность ребер для каждой грани дается обходом против часовой стрелки, если смотреть на тело снаружи. Благодаря этому сохраняется информация о том, с какой стороны граней находится внутренний объем (материал изделия).
Эта же информация может быть использована для описания моделей в системах поверхностного моделирования. В этом случае правило обхода ребер значения не имеет.
Если из табл. 1 убрать сведения о гранях, то оставшаяся информация будет соответствовать системам каркасного моделирования.
Фактически для реальных моделей информация B-Rep более сложна по следующим причинам:
В существующих программных комплексах информация B-Rep организована в виде двусвязных списков, что имеет преимущества перед табличным хранением данных.
Большинство систем твердотельного моделирования основано на структуре данных B-Rep.
Представление кривых и работа с ними
С помощью кривых линий моделируются криволинейные стержни и ребра твердотельных моделей.
Для описания кривых используются уравнения не ниже третьего порядка, потому что они обладают важным свойством: две кривые третьего порядка могут быть соединены друг с другом так, что вторые производные (кривизны) в точке соединения будут равны друг другу. Это означает, что две кривые будут казаться единым целым.
Кривая, отображаемая на экране, в действительности представляет собой набор коротких отрезков. Для их вывода на экран необходимо знать координаты точек кривой, находящихся на равном расстоянии друг от друга. Это называется вычислением кривой.
Эрмитовы кривые
Уравнение Эрмитовой кривой записывается в матричной форме:
.
{x y z} – координаты точки на кривой.
fi(u) – функции сопряжения, которые являются кубическими полиномами относительно параметра u. Параметр u находится в интервале 0 u 1. Значение u=0 соответствует начальной точке кривой, u=1 – конечной.
В прямоугольной матрице по строкам записаны:
Обычно точки в уравнениях кривых записывают более кратко:
.
Вычисление Эрмитовой кривой заключается в подстановке в функции сопряжения параметра u с определенным шагом и нахождении координат точки P.
Кривые Безье
По сравнению с Эрмитовыми кривыми, кривые Безье обеспечивают большее удобство работы (создания алгоритмов и программирования).
Кривая Безье имеет вид:
, где .
Bi,n(u) – функции сопряжения, Pi – образующие точки.
Параметр u полностью аналогичен параметру u Эрмитовых кривых.
Величина n определяет степень кривой Безье.
Рассмотрим для примера кривые Безье третьего порядка. При n=3 получим формулу:
P = (1 – u)3 P0 + 3 u (1 – u)2 P1 + 3 u2 (1 – u) P2 + u3 P3.
На рис. 18 показано графическое представление последней формулы при различных образующих точках при z=0. На рисунке видно, как смещение одной точки изменяет форму кривой.
Существуют выражения для производных кривых Безье. Они необходимы для решения геометрических задач, например, сопряжения кривых, поиска точки пересечения двух кривых и т.д. Производные выводятся дифференцированием выражения P по параметру u.
B-сплайны
Из формулы кривых Безье очевидно, что все задающие точки влияют на форму всей кривой. Это создает следующие неудобства:
В САПр желательно наличие прямо противоположного свойства – локальности изменений.
Локальность изменений предъявляет к кривым следующие требования:
Этим требованиям удовлетворяют B-сплайны. B-сплайн – это составная кривая, состоящая из нескольких разных кривых.
Построение B-сплайна состоит из следующих шагов:
Введем обозначение: k – порядок B-сплайна. В точках сопряжения кривых соблюдается равенство первых k-2 производных. Так, для кубических B-сплайнов равны углы наклона, для B-сплайнов четвертой степени – кривизны и т.д.
Существуют алгоритмы деления и объединения B-сплайнов, причем в результате этих действий порядок k не изменяется.
Интерполяционные кривые
К интерполяционным кривым относятся кривые, уравнения которых определяются по заданным точкам, лежащим на самой кривой.
Рассмотрим интерполяцию на примере Эрмитовых кривых.
Исходными данными являются точки P0, P1, P2, … , Pn (рис. 19).
Для определения кривых необходимо найти углы наклона P0’, P1’, P2’, … , Pn’ в каждой точке.
Таким образом, задача имеет n+1 неизвестное. Задача будет решена, если удастся записать n+1 независимое уравнение.
Предполагается, что общая кривая имеет свойство равенства кривизн в точках сопряжения, то есть (первый индекс обозначает номер точки, второй – номер кривой). Это дает n-1 уравнение.
Недостающие два условия определяют одним из двух следующих способов.
После нахождения углов наклона P0’, P1’, P2’, … , Pn’ Эрмитовы кривые полностью определены.
Представление поверхностей и работа с ними
Для представления поверхностей пользуются такими же обозначениями, как и для представления кривых:
Рассмотрим некоторые способы аналитического описания поверхностей.
Билинейная поверхность
Билинейная поверхность строится по четырем заданным угловым точкам (рис. 20):
.
Вдоль каждого параметра (u, v) вертикальная координата изменяется линейно, и линии сетки на поверхности – прямые. Однако в целом билинейная поверхность будет являться плоскостью только тогда, когда все четыре угловые точки будут находиться в одной плоскости.
Бикубический лоскут
Бикубический лоскут определяется уравнением
.
Можно заменить исходные алгебраические коэффициенты aij параметрами, имеющими для поверхности геометрический смысл. Для этого необходимо составить 16 уравнений (по количеству коэффициентов aij), которые задают геометрические параметры поверхности в характерных точках. Такими уравнениями являются следующие условия.
Система уравнений решается относительно aij, выражения aij подставляются в формулу P, в результате чего уравнение бикубического лоскута принимает вид:
.
Функции сопряжения fi полностью совпадают с функциями сопряжения Эрмитовых кривых. Таким образом, бикубический лоскут – это, во-первых, расширение Эрмитовой кривой на одно измерение, а во-вторых, описание семейства Эрмитовых кривых, соответствующих различным значениям v.
Поверхность Безье
Можно расширить на одно измерение концепцию кривой Безье, в результате чего получится поверхность Безье:
,
где Pij – задающие точки, находящиеся в вершинах задающего многогранника;
Bin(u), Bjm(v) – функции сопряжения кривых Безье.
При моделировании сложных поверхностей получаются уравнения Безье высоких степеней, что создает проблемы вычислительного характера. Поэтому обычно используют уравнения Безье степени 3 по u и v на небольших участках поверхности, которые соединяют друг с другом.
Аналогично кривым Безье, существуют формулы производных для поверхностей. Расчет по этим формулам требуется, например, при решении задач тонирования, когда требуется вычислить угол падения света, и в некоторых других случаях.
B-сплайновые и интерполяционные поверхности
Подобно тому, как от уравнения кривой Безье можно перейти к уравнению поверхности Безье, от уравнения B-сплайна и интерполяционной кривой можно перейти к уравнениям B-сплайновой поверхности и интерполяционной поверхности.
С помощью B-сплайновых и интерполяционных поверхностей представляются поверхности, полученные методами охвата по шаблону, скиннинга, скругления, плавного сопряжения и моделирования границ.
МАТЕМАТИЧЕСКОЕ ОБЕСПЕЧЕНИЕ МОДУЛЯ CAE
Созданная геометрическая модель должна быть проанализирована в условиях эксплуатации, когда она подвергается силовым и прочим воздействиям. Моделирование условий эксплуатации происходит по алгоритмам функционирования объектов при каких-либо внешних воздействиях.
Описание функционирования объектов производится с помощью алгебраических, дифференциальных и других уравнений. Системы таких уравнений решаются численными методами.
Алгебраические уравнения
Многие расчетные задачи приводятся к системе n линейных уравнений с n неизвестными вида [A] {X} = {B}, в которой матрица [A] является симметричной и ленточной. Ленточной называется такая матрица, у которой все элементы вне полосы, расположенной вдоль главной диагонали, равны нулю.
Задача [A] {X} = {B} разрешима, если матрица [A] является неособенной (регулярной), т.е. если ее определитель не равен нулю.
В противном случае матрица [A] называется особенной (сингулярной).
Числом обусловленности матрицы [A] называется величина, равная произведению определителей
C(A)= ║A║║A-1║.
Если матрица [A] становится сингулярной, то C(A) +.
Если C(A) +, то матрица [A] является плохо обусловленной.
Численное решение системы [A] {X} = {B} с плохо обусловленными матрицами затруднено тем, что малые возмущения исходных данных ведут к большим возмущениям решения.
При решении систем с такими матрицами метод Гаусса (и ему подобные методы) неудобен.
В таких случаях возможно применение итерационных методов. Последовательно вычисляются приближения {X}0, {X}1, … , {X}k к неизвестному решению по формулам типа {X}i=f({X}i-1). На каждом шаге проверяется степень совпадения левой ([A] {X}i) и правой ({B}) частей уравнения.
В общем случае не все итерационные методы обладают свойством сходимости при произвольных исходных данных.
Однако итерационные методы обеспечивают сходимость в некоторой окрестности решения, которая может быть очень малой. Поэтому скорость решения (и даже сама возможность получения решения) зависит от удачного или неудачного выбора начального приближения {X}0.
Дифференциальные уравнения с начальными условиями
Пусть требуется решить задачу , если известно, что искомая функция y = f(x) проходит через точку (x0, y0).
В общем случае используется численный метод, состоящий в замене производной конечными разностями. Эту замену можно провести различными способами. Простейший из них приводит к явной формуле Эйлера (см. курс лекций и методические указания к практическим занятиям по дисциплине «Математическое моделирование»).
На практике применяют более сложные формулы, например, формулы Рунге-Кутта.
Проекционные методы
Рассмотрим уравнения, описывающие напряженно-деформированное состояние конструкции:
L[y(x)] = f(x),
где L[y(x)] – дифференциальный оператор над функцией y(x).
Пример.
Предположим ситуацию, когда разрабатывается новый конечный элемент в виде стержня с шарнирами по концам, загруженного в середине пролета сосредоточенной силой P. Требуется найти уравнение изогнутой оси стержня в зависимости от параметров L (пролета), P (нагрузки), EI (жесткости).
Дифференциальное уравнение изгиба стержня имеет вид M=EIy’’, где y=y(x) – функция, описывающая изогнутую ось стержня.
Уравнение равновесия произвольного сечения стержня:
|
M = (P/2) x |
для x = 0 … L/2; |
|
M = (P/2) (L-x) |
для x = L/2 … L. |
Приравняв величины изгибающих моментов, получим:
|
EI y’’ = (P/2) x |
для x = 0 … L/2; |
|
EI y’’ = (P/2) (L-x) |
для x = L/2 … L. |
В полученных формулах: L[y(x)] = EI y’’; функция f(x) представлена правыми частями.
Решение y(x) удобно искать в виде линейной комбинации независимых функций (x), каждая из которых удовлетворяет граничным условиям:
y(x) = Ai i(x).
Обычно в качестве функций (x) берут тригонометрические функции. После выбора функций (x) решение задачи сводится к определению неизвестных коэффициентов Ai. Точность решения будет зависеть от правильного выбора функций (x) и от их количества в сумме Ai i(x).
Например, решение в виде проиллюстрировано на рис. 21. Можно подобрать такие слагаемые , что изогнутая ось стержня с хорошей точностью будет соответствовать любой нагрузке.
Все рассмотренные ниже методы пригодны для решения различных дифференциальных уравнений в одно-, двух- и трехмерном пространстве.
Получаемые решения будем сравнивать с точными значениями прогибов стержня:
Введем понятие «невязки»: (x) = L[y(x)] - f(x). Невязка используется для реализации некоторых методов.
Метод коллокации
Метод коллокации заключается в приравнивании нулю невязки в заданных точках. Количество точек (n) должно совпадать с количеством учитываемых функций (x).
Рассмотрим вариант n=1.
Зададимся искомой функцией в виде и найдем коэффициент A1.
Приравняем нулю невязку в середине пролета:
, или .
После подстановки в последнюю формулу x=L/2 и решения уравнения относительно A1 получим:
;
.
Вычислим прогибы в четверти и середине пролета по последней формуле.
y(x=L/4) = 0,85975 y0. Погрешность (относительно 0,6875 y0) составляет 25,1%.
y(x=L/2) = 1,2159 y0. Погрешность составляет 21,6%.
Рассмотрим вариант n=2.
(взяты симметричные функции, удовлетворяющие граничным условиям, при минимальных целых коэффициентах перед x/L).
Выражение для невязки принимает вид: .
Для нахождения коэффициентов A1 и A3 подставим в это уравнение величины x=L/4 и x=L/2. Получим систему двух линейных алгебраических уравнений относительно A1 и A3, решение которой дает:
A1 = 1,0378 y0;
A3 = 0,0198 y0;
.
Вычислим прогибы в четверти и середине пролета.
y(x=L/4) = 0,7198 y0. Погрешность составляет 4,7%.
y(x=L/2) = 1,0576 y0. Погрешность составляет 5,8%.
Рассмотрим, заведомо худший вариант при n=2. Он иллюстрирует, каким образом неудачно выбранные функции влияют на точность расчетов.
Возьмем функцию . Второе слагаемое не удовлетворяет граничным условиям: при x=0 и x=L оно дает ненулевые прогибы.
Выражение для невязки: .
Подстановка в уравнение x=L/4 и x=L/2 приводит к решению:
A1 = 0,85973 y0;
A2 = 0,08903 y0;
.
Прогибы на опоре, в четверти и середине пролета:
y(x=0) = 0,08903 y0, что не соответствует условиям задачи;
y(x=L/4) = 0,60792 y0 (погрешность составляет 13,1%);
y(x=L/2) = 0,94876 y0 (погрешность составляет 5,1%).
Метод наименьших квадратов
Составляется аналитическое выражение . Для вычисления коэффициентов Ai решается система уравнений . Размерность системы (nn) равна количеству функций в сумме Ai i(x).
Рассмотрим вариант n=1.
Функция прогибов: .
Для получения функции Q необходимо проинтегрировать выражение [(x)]2 на интервале от 0 до L. При этом необходимо учесть, что выражение для невязки различно при x L/2 и x L/2.
.
После возведения подынтегральных выражений в квадрат, интегрирования и алгебраических упрощений получим:
.
Запишем уравнение : и решим его относительно A1: .
Искомая функция: .
Прогибы в четверти и середине пролета:
y(x=L/4) = 0,69687 y0 (погрешность составляет 1,36%);
y(x=L/2) = 0,98553 y0 (погрешность составляет 1,45%).
Метод моментов
Записываются уравнения ортогональности невязки (x) к x0, x1, …, xn-1, где n – количество слагаемых в сумме Ai i(x). Уравнение ортогональности в общем виде можно записать как , где k = 0, 1, …, n-1.
Рассмотрим вариант n=1.
Функция прогибов: .
Уравнение ортогональности:
.
После интегрирования и алгебраических упрощений получим:
.
Решение уравнения относительно A1: .
Искомая функция: .
Прогибы в четверти и середине пролета:
y(x=L/4) = 0,67542 y0 (погрешность составляет 1,78%);
y(x=L/2) = 0,95493 y0 (погрешность составляет 4,51%).
Метод Галёркина
Записываются уравнения ортогональности невязки (x) к первым n функциям i(x): .
Рассмотрим вариант n=1.
Функция прогибов: , функция .
Уравнение ортогональности:
.
После алгебраических упрощений и интегрирования получим решение, совпадающее с решением по методу наименьших квадратов.
Рассмотрим вариант n=2.
Функция прогибов: .
Функция , функция .
Можно составить два уравнения ортогональности:
После алгебраических упрощений, интегрирования и решения уравнений относительно A1, A3 получим: , .
Искомая функция: .
Прогибы в четверти и середине пролета:
y(x=L/4) = 0,68827 y0 (погрешность составляет 0,11%);
y(x=L/2) = 0,9977 y0 (погрешность составляет 0,23%).
Метод Галеркина часто оказывается точнее других методов, его сходимость лучше изучена, и он служит основой метода конечных элементов.
Метод Рэлея-Ритца
Записывается выражение для полной (внутренней плюс потенциальной) энергии T системы. Известно, что искомая функция y(x) минимизирует T. Это дает возможность найти неизвестные коэффициенты Ai путем решения системы уравнений .
Полная энергия стержня, изгибаемого сосредоточенной силой, равна: .
Первое слагаемое выражает накопленную энергию изгиба, второе – работу внешней силы P на перемещении yc в точке ее приложения (середине стержня).
Рассмотрим вариант n=1.
Функция прогибов: . Прогиб в середине стержня (при x=L/2) равен A1. С учетом этого полная энергия выразится уравнением:
.
После дифференцирования функции y(x), алгебраических преобразований, интегрирования, решения относительно A1 уравнения получим решение, совпадающее с решением по методу наименьших квадратов и методу Галеркина.
Рассмотрим вариант n=2.
Функция прогибов: . Подставив в функцию прогибов значение x=L/2, получим, что прогиб в середине стержня равен A1A3. С учетом этого полная энергия выразится уравнением:
.
После дифференцирования функции y(x), алгебраических преобразований, интегрирования, решения относительно A1 и A3 системы уравнений , получим решение, совпадающее с решением по методу Галеркина.
Метод конечных разностей
Этот метод состоит в замене производных приближенными выражениями, в которые входят значения функций в некоторых точках.
На рис. 22 показана геометрическая интерпретация приближенного выражения первой производной некоторой функции, и даны обозначения, используемые в приближенных выражениях производных более высоких порядков.
Первую производную в точке i можно приближенно найти как тангенс угла наклона отрезка, соединяющего соседние точки:
.
Запишем аналогичные выражения для первых производных на участках a, b, примыкающих к точке i: , .
Так как , то , что приводит к формуле
.
Записав первые производные , , и приняв во внимание то, что , получим формулу
.
Записав вторые производные , , и приняв во внимание то, что , получим формулу
.
Разобьем стержень на 4 участка (рис. 23, а) и будем искать прогибы y1, y2, y3 в точках 1, 2, 3. В силу симметричности задачи можно рассматривать только половину стержня. Запишем уравнение EI y’’ = (P/2)x для точек 2, 3 с использованием приближенной формулы второй производной:
В двух уравнениях имеются три неизвестных. Третье уравнение, необходимое для решения задачи, можно записать, используя граничные условия. В данном случае таким уравнением является очевидное условие y1=0.
Решением уравнений является:
. Погрешность составляет 9,1%.
. Погрешность составляет 12,5%.
Попробуем получить более точное решение и разобьем стержень на 8 участков (рис. 23, б). Система уравнений принимает вид:
Система имеет решение относительно y1 … y5:
y2 = 0,375 y0.
y3 = 0,703125 y0. Погрешность равна 2,3%.
y4 = 0,9375 y0.
y5 = 1,03125 y0. Погрешность равна 3,1%.
Метод конечных разностей прост в реализации, легко адаптируется к различным дифференциальным уравнениям и уравнениям в частных производных.
Использование метода конечных разностей ограничено трудностью наложения регулярной квадратной сетки на области произвольной формы.
Метод конечных элементов (МКЭ)
МКЭ применяется для решения самых разнообразных задач: в строительной механике, для решения задач теплопереноса, распределения электростатического потенциала, механики жидкостей, вибрационного анализа и т.д.
При переходе от геометрической модели к модели, пригодной для расчетов методом КЭ, проектировщику приходится принимать одно из наиболее важных решений – разбиение модели на конечные элементы и выбор подходящих элементов с нужным количеством узлов из библиотеки доступных элементов.
Разбиение модели на конечные элементы называется построением сетки.
При моделировании стержневых систем задача упрощается тем, что независимо от длины конечных элементов будет получено точное решение. Следовательно, конечноэлементная стержневая модель должна содержать минимальное количество узлов и элементов. Узлы необходимо вводить, например, в точки изменения направления продольной оси элемента, в точки примыкания других элементов, в места изменения жесткостных характеристик и т.д.
При моделировании поверхностных или твердотельных систем построение сетки относится к наиболее сложному этапу моделирования. В указанных системах точность решения зависит от величины конечных элементов. Проектировщику приходится учитывать некоторые важные моменты, например:
Многие САПр поддерживают методы автоматического построения сетки. Построение сетки подразумевает определение положения узлов и элементов на основании минимального объема данных, вводимых проектировщиком.
Рассмотрим существующие методы построения сетки.
Создание узлов методом Кавендиша
Работа начинается с задания граничных узлов вручную (рис. 24). Далее строится квадратная сетка с заданным размером ячеек. В каждом квадрате сетки случайным образом создается один узел. Если узел попадает внутрь объекта, и расстояние от него до граничных и ранее созданных узлов превышает размер ячейки, этот узел считается принятым. В противном случае в том же квадрате случайным образом создается другой узел, который проходит ту же проверку. Если за фиксированное количество попыток (например, 5) принять узел не удается, происходит переход к следующему квадрату.
Этот метод может быть расширен до трех измерений простым переходом от плоской квадратной сетки к пространственной кубической.
Создание узлов методом Шимады
В центре заданных граничных узлов помещаются опорные шары, на которые натягивается граничная поверхность (рис. 25). Полученный объем заполняется шарами так, как будто на них действует сила тяжести, или сила, направленная от центра к границам, либо используется какой-либо другой алгоритм заполнения объема. Размер шаров подбирается по заданному максимальному размеру конечных элементов. Центры шаров становятся узлами.
Построение элементов
Наиболее популярным методом соединения узлов является метод триангуляции Делоне. Данный метод ориентирован на то, чтобы по возможности избегать формирования узких треугольников.
Триангуляция Делоне реализуется несколькими способами. В том числе, она может производиться непосредственно на наборе узлов с использованием алгоритма двумерной триангуляции Ватсона. Согласно этому алгоритму, три узла, не лежащие на одной прямой, объединяются в треугольник, если окружность, проведенная через эти три узла, не охватывает других узлов.
Рассмотрим фрагмент построения элементов по узлам, полученным методом Шимады (рис. 26). Предположим, узлы отсортированы по возрастанию вертикальной координаты и соответственно этому пронумерованы. Происходит перебор троек узлов из полного их множества, например, 1-2-3, 1-2-4, 1-2-5 … и т.д. Через каждую тройку узлов проводится окружность (для объемных моделей – сфера), и проверяется попадание в нее других узлов. На рис. 26 видно, что окружность, проведенная по узлам 1-2-3, охватила другие узлы. Тем самым удалось избежать формирования узкого треугольника 1-2-3. Тупоугольный треугольник 1-2-7 также не будет создан. Первый треугольный элемент будет создан для сочетания узлов 1-3-7.
Топологическое разбиение и геометрическое разбиение
Это два аналогичных метода, суть которых заключается в том, что область разбиения последовательно делится на все более мелкие элементы до достижения заданного размера конечных элементов.
Решеточные методы
Решеточные методы основаны на том, что на область разбиения накладывается сетка. Если ячейка сетки полностью попадает в область разбиения, то на месте ячейки формируется конечный элемент (рис. 27). Ячейки сетки вдоль границ преобразуются в треугольные конечные элементы методом триангуляции. Разновидности решеточных методов отличаются друг от друга главным образом методами создания граничных элементов, а также тем, что сетка может быть прямоугольной и треугольной.
К решеточным методам относится представление октантного дерева при декомпозиционном моделировании. На завершающем этапе граничные октанты разбиваются на тетраэдры для более точного моделирования границ.
Отображаемые элементы
На рис. 28 показано разбиение на конечные элементы области произвольной формы. Исходная область отображается на регуляризованный прямоугольник (или параллелепипед в пространстве). На регуляризованную область наносится сетка в соответствии с заданными размерами конечных элементов. Сетка отображается на исходную область с образованием конечных элементов. Конечные элементы получаются более сложными и требуют больших вычислительных ресурсов.
Существует много методов отображения, среди которых самыми типичными являются трансфинитное и изопараметрическое отображение.
При трансфинитном отображении обрабатываются аналитические уравнения границ исходной области. В результате сетка конечных элементов строится без искажения границ.
При изопараметрическом отображении соответствие между исходной и регуляризованной областями обеспечивается в нескольких граничных точках. Исходные границы заменяются сплайнами, проходящими через указанные точки.
Оптимизация
В принципе весь процесс проектирования можно считать оптимизацией в широком смысле, поскольку проектировщик стремится получить проектное решение, которое можно реализовать с минимальными затратами.
Обычно (в более узком смысле) под оптимизацией понимается выбор оптимальных размеров конструкции.
Оптимизация требует параметризации конструкции, дающей возможность рассматривать семейство похожих конструкций, изменяя значения параметров. Например, параметрами фермы, изображенной на рис. 29, являются пролет L, высота h, размер панели верхнего пояса a, тип решетки, нагрузка q, поперечные сечения элементов, характеристики материала. Различные наборы значений параметров будут давать различные фермы.
В зависимости от ситуации некоторые параметры не могут варьироваться. Обычно величины L, a, q в расчетах участвуют как константы. Другие параметры являются оптимизируемыми и называются переменными оптимизации.
Показателем качества фермы может быть отношение . Показатель качества, вычисляемый по переменным оптимизации, называется целевой функцией.
Обозначим переменные оптимизации символом {X} (n-мерный вектор, компонентами которого являются переменные оптимизации), а целевую функцию символом F({X}).
Очень редко показатель качества может быть выражен одной целевой функцией. Чаще всего приходится выбирать между разными показателями или строить объединенный показатель с какими-либо весовыми коэффициентами. Этот процесс называется построением сложной целевой функции.
В большинстве случаев переменные оптимизации имеют ограниченную область определения. Например, расчетное сопротивление материала не выходит за пределы известного интервала. Ограничение, задающее границы области определения переменной оптимизации, называется ограничением области.
Ограничение, полученное из явного рассмотрения функционального требования к конструкции, называется функциональным ограничением. Например, напряжения в элементах фермы не должны превышать расчетного сопротивления материала.
Таким образом, простая задача оптимизации может быть записана следующим образом:
найти {X}*, такой что F({X}*) = min F({X}) или F({X}*) = max F({X})
при условии, что
xi1 xi xi2 (ограничения области на i-ю переменную);
Gk({X}*) 0 (функциональные ограничения).
Методы поиска
Методы поиска экстремума целевой функции могут быть сгруппированы в три больших класса (рис. 30). На рис. 30 не указаны простейшие методы поиска (по одной переменной оптимизации).
Рассмотрим наиболее употребительные методы.
Методы, основанные на градиентах, используют первые производные. Методы, основанные на матрице Гессе, используют матрицу, составленную из значений самой функции, первых производных и вторых частных производных . Матрица Гессе в общем случае требует достаточно ресурсоемких вычислений.
Прямые методы используют сведения о локальном поведении функции и сходятся к локальному экстремуму. С их помощью трудно анализировать все пространство поиска.
Методы направленного случайного поиска просматривают все пространство оптимизации, поэтому с их помощью можно найти глобальный экстремум целевой функции. Они все шире используются во многих приложениях.
Если оптимизацию строительной конструкции привести к комбинаторной оптимизации, то можно найти решение, близкое к оптимальному. Для фермы, изображенной на рис. 29, тип решетки, поперечные сечения элементов, расчетное сопротивление материала являются дискретными переменными. Если задать изменение высоты в каких-либо пределах с определенным шагом, то получим задачу комбинаторной оптимизации.
Предположим, что ищется минимум целевой функции и рассмотрим алгоритм комбинаторной оптимизации.
1) Назначается начальное значение вектора {X}, которое принимается за искомое решение {X}*. Вычисляется целевая функция F({X}*).
2) Каждая переменная оптимизации изменяется на величину kS. S – шаг изменения переменной, k – величина, случайным образом принимающая значения -1, 0, 1. Для нового вектора {X} вычисляется целевая функция F({X}).
3) Если = F({X}) – F({X}*) 0, то {X}* = {X} (за решение задачи принимается текущий вектор {X}).
4) Если = F({X}) – F({X}*) 0, то {X}* = {X} с некоторой вероятностью p, зависящей от .
5) Переход к п. 2.
Выход из алгоритма и завершение поиска происходит в п. 4, когда за заданное количество попыток вектор {X}* не изменяется.
Метод модельной закалки разрешает отклонения от минимума, но дает возможность попасть в другой минимум, более глубокий.
Применение метода к конкретной задаче требует аккуратности в задании исходных данных (начального вектора {X} и шага S каждой переменной), а также зависимости p().
Рассмотрим получение начального состояния системы.
Назначаются M начальных значений вектора {X}. Предположим, целевая функция вычисляется по N переменным оптимизации. Каждая переменная кодируется двоичной строкой. Двоичное представление переменной определяет ген информации. Набор генов, соответствующий N переменным, определяет хромосому. Набор хромосом – это популяция двоичных строк, которая будет развиваться по принципам естественного отбора. Каждой хромосоме ставится в соответствие степень жизнеспособности – величина целевой функции, вычисленной по значениям переменных оптимизации. Начальное состояние системы показано на рис. 31.
Для начального кодирования переменных широко используется метод целочисленного представления. Предположим, переменная принимает значения от –/2 до +/2, которые необходимо представлять с точностью до двух знаков после запятой. Получим область определения переменной [–1,57; 1,57]. Умножением на 100 значения переменной отображаются на целые числа от –157 до +157. Для их представления необходимы 8 битов плюс 1 бит на знак числа. На рис. 31 ген11 равен значению +157.
Рассмотрим механизмы воспроизводства и отбора.
1) Определение родителей. Этот этап представлен в табл. 2 на примере четырех хромосом и максимизации целевой функции.
Таблица 2 – Определение родителей
|
Хромосома |
F({X}) |
Вероятность pi = Fi / Fm |
Количество потомков |
|
|
X1 |
165 |
0,574 |
0,574 |
|
|
X2 |
570 |
1,983 |
1 + 0,983 |
|
|
X3 |
65 |
0,226 |
0,226 |
|
|
X4 |
350 |
1,217 |
1 + 0,217 |
|
|
Сумма Fi |
1150 |
4 |
4 |
|
|
Среднее Fm = Fi / 4 |
287,5 |
Для каждой хромосомы вычислена целевая функция и вероятность участия в образовании новых хромосом. Если вероятность больше единицы, хромосома (в данном случае X2 и X4) участвует в образовании одной новой хромосомы (потомка), и из ее вероятности единица вычитается. Оставшееся количество потомков распределяется между имеющимися хромосомами (в данном случае X1 и X2) по убыванию вероятностей.
2) Воспроизводство. Из хромосом-родителей (в предыдущем примере – X1, X2, X2, X4) случайным образом выбираются пары хромосом. Каждая пара разрезается на части в случайных местах, происходит обмен частями, в результате чего образуются новые хромосомы (рис. 32). Этот процесс называется кроссовером (cross over). Кроссовер выполняется с некоторой вероятностью. Если кроссовера не происходит, то родительские хромосомы просто дублируют сами себя.
3) Мутация. Мутация применяется индивидуально к каждому потомку после кроссовера. Оператор мутации с малой вероятностью (например, 0,001) изменяет каждый ген путем замены какого-либо 0 на 1 или наоборот.
4) Отбор. После пополнения популяции потомками они должны быть проверены на соответствие ограничениям типа G({X}) 0. Потомки, не удовлетворяющие ограничениям, из популяции удаляются. Кроме того, из популяции могут удаляться хромосомы с низкой степенью жизнеспособности.
5) Переход к п.1.
При правильной реализации генетического алгоритма популяция развивается с увеличением степени жизнеспособности хромосом, и целевая функция стремится к глобальному оптимуму.
Структурная оптимизация
Структурной оптимизацией называется автоматический синтез компонентов модели.
В зависимости от того, какие свойства компонентов оптимизируются, бывает оптимизация размеров, формы или топологии.
Оптимизация размеров
Оптимизация размеров состоит в изменении размеров конструкции при сохранении ее формы и топологии. Для фермы с фиксированным пролетом (рис. 29) оптимизация размеров заключается в определении оптимальных поперечных сечений элементов.
К оптимизации размеров можно отнести также выбор оптимального материала.
Оптимизация формы
Оптимизация формы подразумевает сохранение неизменной топологии при изменении формы. На рис. 33 показана ферма, которая получена из фермы, изображенной на рис. 29, путем изменения формы с сохранением топологии.
В этом случае в качестве переменных оптимизации можно использовать координаты узлов фермы (не обязательно всех узлов), либо важнейшие размеры фермы.
Оптимизация топологии
Оптимизация топологии связана с изменением внутренней структуры модели, с созданием и удалением компонентов модели. На рис. 34 показана ферма, которая отличается топологией от фермы, изображенной на рис. 29.
Рассмотрим существующие подходы к оптимизации стержневых систем. Общим в существующих подходах является то, что оптимизация требует многократного выполнения расчетов методом конечных элементов.
Пространство делится на одинаковые прямоугольные ячейки. Внутри каждой ячейки имеется пора размерами xi yi zi, где i – номер ячейки. Размеры пор совместно с углом ориентации ячеек рассматриваются как переменные оптимизации.
Алгоритм оптимизации увеличивает размеры полостей в тех областях, где материал нагружен недостаточно сильно, и наоборот. Ориентация ячеек подбирается таким образом, чтобы жесткость системы была максимальной.
Литература