C(p(t)) = 0.

Целевая конфигурация задается указанием желательных значений для некоторых координат в конфигурационном пространстве: для прямой задачи - требуемых значений переменных управления кинематическими парами, для обратной задачи - требуемых значений переменных, задающих положение нужных звеньев в пространстве. В целях получения натурального решения логично задать целевые значения и для всех остальных переменных - равные их начальным значениям. Тем самым будет гарантировано нахождение минимального решения (самого короткого движения механизма). В этих условиях положение целевой конфигурации в конфигурационном пространстве задается точкой p*. Заметим, что в общем случае C(p*) 0. Но даже если целевая конфигурация удовлетворяет всем кинематическим связям, для решения задачи кинематики требуется найти траекторию движения механизма p(t).

Дифференциальное уравнение движения

Наиболее натуральным способом моделирования движения механизма является представление его в виде решения дифференциального уравнения. В каждой точке траектории осуществляется линеаризация системы уравнений C(p(t)):

C(p(t + dt) = JC(p)pdt,

где JC(p) - матрица Якоби (частных производных) системы уравнений C, вычисленная в точке p. В предположении C(p(t)) 0 для всех t получаем JC(p)p = 0, то есть . Предположим, M - ортонормированный базис Ker JC(p). Тогда линейный оператор проектирования на это подпространство будет иметь вид MMT, а дифференциальное уравнение кинематики запишется как

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

Натуральный градиент уравнения

Для решения дифференциального уравнения, описанного в предыдущем параграфе, необходимо научиться дифференцировать уравнения кинематических пар в произвольной точке p. Эти уравнения связывают координаты ключевых точек, прямых и плоскостей звеньев с помощью геометрических и инженерных ограничений. Геометрические координаты объектов вычисляются по их начальным координатам путем применения трансформаций, задающих положение жестких множеств (звеньев) в трехмерном пространстве. Множество трансформаций образуют группу (с операцией композиции), называемую специальной евклидовой группой SE(E3). На этой группе определена метрика (в виде евклидовой нормы разности матриц трансформаций), а значит, она является группой Ли. Таким образом, нам удобней заменить введенную в предыдущем параграфе функцию уравнения Ci: Rn R на сложную функцию Ci(Dk(p),Dl(p)), где Dj: Rn х SE(E3), а Ci: SE(E3) SE(E3) х R. Таким образом, вектор-функция D: Rn х (SE(E3))k описывает положения всех k звеньев с помощью трансформаций, а каждое уравнение задает связь между двумя звеньями, используя эти трансформации. Натуральным градиентом уравнения Ci: SE(E3) х SE(E3) х R назовем вектор С(T1, T2) = (a, b, c, d) х R12, где

Функция exp задает экспоненциальную параметризацию однопараметрической подгруппы группы Ли SE(E3). Натуральный градиент позволяет вычислить частные производные функций Ci по переменным конфигурационного пространства с использованием операций алгебры Ли.

Алгоритмы численного решения дифференциальных уравнений

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

? метод Эйлера;

? классический метод Рунге-Кутта;

? метод Дорманда-Принса (DOPRI5).

Все эти методы являются частными случаями общей схемы Рунге-Кутта для решения дифференциального уравнения с начальным условием:

где y: R Rn - искомая, а f: Rn+1 Rn - известная функции. Общая схема Рунге-Кутта определяется так:

Здесь s задает порядок схемы, h - величину шага, а константы aij определяют конкретный метод. При s = 1 используется метод Эйлера, решающий дифференциальное уравнение с точностью третьего порядка (от величины шага), при s = 2 работает классический метод Рунге-Кутта (обладающий точностью четвертого порядка), а при s = 7 применяется метод Дорманда-Принса с точностью пятого порядка. На практике предпочтительней применять последний метод, несмотря на то, что он требует вычисления на каждом шаге большего количества значений.

Планирование движения

Задача планирования движения задается описанием:

? окружения;

? устройства (манипулятора, робота);

? начальной и целевой конфигураций устройства.

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

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

Рис. 30

Таким образом, для построения дорожной карты необходимо уметь выполнять три действия:

? проверить заданную конфигурацию на отсутствие пересечений частей устройства и препятствий;

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

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

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

Вопросы для самоконтроля

1. Дайте определение прямой и обратной задачам кинематики.

2. Опишите основные кинематические пары.

3. Как моделируются механизмы в терминах задач удовлетворения ограничениям?

4. Что такое конфигурационное пространство механизма?

5. Каким дифференциальным уравнением описывается движение механизма?

6. Что такое натуральный градиент уравнения кинематической пары, и для чего он нужен?

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

8. Как осуществляется планирование движения с помощью дорожной карты?

Дополнительная литература

Хорошее введение в предметную область кинематики механизмов можно получить из учебников по технической механике, например [12]. Представление о методах параметризации группы SE(E3) и дифференцирования ее однопараметрических подгрупп можно получить из работы [38]. Использование групп Ли для расчета движения механизма в конфигурационном пространстве описано в [17]. Обзор известных методов планирования движения дан в [31].



Лекция 11. Инженерный анализ динамики

Задача анализа динамики механизмов

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

Движение абсолютно твердого тела в трехмерном пространстве

Как известно, движение частицы с массой m в поле действия сил {fi}i=1,…,n описывается вторым законом Ньютона, постулирующим равенство изменения импульса частицы (определяемого произведением ее массы и вектора скорости) действующим на нее силам. В предположении неизменной массы частицы, второй закон Ньютона выглядит так:

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

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

где общая масса тела M определяется интегралом плотности по объему:

Обозначив скорость центра масс vc, получаем

Уравнение вращения тела вокруг его центра масс получается интегрированием уравнения изменения момента импульса для точечных масс:

Где - угловая скорость, а I - тензор инерции тела:



Таким образом, из уравнения изменения момента импульса следует

где ri - точка приложения к телу силы fi

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

Моделирование контакта тел

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

? пружинно-амортизаторные пары (метод штрафов);

? двусторонние контактные силы;

? односторонние контактные силы;

? повторный импульс.

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

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

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



Альтернативный подход: уравнения Лагранжа

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

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

Пусть в определенные моменты времени t = t1 и t = t2 система занимает определенные положения, характеризуемые двумя наборами координат q(1) и q(2). Тогда между этими положениями система движется таким образом, чтобы интеграл имел наименьшее возможное значение.

Функция L называется функцией Лагранжа данной системы, а вышеприведенный интеграл - действием. Нетрудно вывести, что минимизация этого интеграла эквивалентна решению системы дифференциальных уравнений, называемых уравнениями Лагранжа:

или в векторном виде:

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

Функция Лагранжа твердого тела равна разности его кинетической и потенциальной энергии:

Заметим, что



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

где J - матрица Якоби системы уравнений, задающей кинематические связи в терминах обобщенных координат, а - вектор неизвестных множителей Лагранжа. Получаемая система содержит меньше уравнений, чем при использовании уравнений Ньютона-Эйлера. Однако для ее решения мы должны также учесть собственно уравнения кинематических связей, то есть решить систему дифференциально-алгебраических уравнений, что является более сложной вычислительной задачей, чем решение системы обыкновенных дифференциальных уравнений.

Методы определения столкновений

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

Заметим, что методы определения столкновений используются в разных приложениях САПР и в общем случае отвечают на один из трех вопросов:

? есть ли столкновение в системе;

? каковы точки контакта и нормали в них;

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

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

Алгоритмы широкой фазы

Все известные алгоритмы широкой фазы в той или иной степени основаны на использовании деревьев ограничивающих объемов (Bounding Volume Trees, BV-trees), представляющих собой аппроксимацию объемного тела набором непересекающихся канонических объемов (таких как параллелепипеды или сферы).

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

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

AABB-деревья окружают каждое тело параллелепипедом со сторонами, параллельными координатным осям (AABB - Axis Aligned Bounded Boxes). Для поиска пересечения параллелепипедов используются их проекции на координатные оси (представляющие собой интервалы). Для каждой оси эти интервалы сортируются, что позволяет обнаружить пересечения. При обнаружении пересечения между двумя параллелепипедами по всем трем координатным осям алгоритм может разделить параллелепипед на четыре части и продолжить сравнение. Обычно используют два разных типа AABB-деревьев: параллелепипеды фиксированных и динамических размеров. Первые окружают тело вокруг его центра масс с учетом всех возможных вращений. Такие параллелепипеды могут быть построены заранее, а затем использованы на каждом шаге алгоритма динамической симуляции движения. Динамический размер подразумевает построение наименьшего параллелепипеда для данной ориентации тела, что требует дополнительных вычислений на каждом шаге симуляции.

OBB-деревья отличаются от AABB тем, что параллелепипеды допускают произвольную ориентацию (OBB - Oriented Bounded Boxes), тем самым каждое тело окружается наименьшим параллелепипедом, который затем вращается вместе с телом. Проверка на пересечение для OBB-деревьев выполняется сложнее, чем для AABB.

K-DOPs (от DOP - Discrete Oriented Polytopes) определяют семейства выпуклых многогранников, грани которых имеют только конечное число (K) ориентаций (нормалей). Например, AABB-деревья являются 6-DOP. В данном подходе подразумевается статическое построение многогранников до начала динамической симуляции движения, а на каждом шаге проверки столкновений вычисляется пересечение заранее построенных K-DOPs вокруг текущих центров масс тел.

Заметаемые сферой объемы (SSV - Swept Sphere Volumes) определяют семейство таких канонических объемов, как сфера, конечный цилиндр со скругленными окончаниями, параллелепипед со скругленными гранями. Сложность их обработки сопоставима со сложностью обработки OBB-деревьев, но данное представление во многих случаях позволяет провести более аккуратную проверку на столкновение.

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

Алгоритмы фазы сужения

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

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

Алгоритм Лина-Канни (его авторы - Ming C. Lin и John F. Canny) основан на свойстве когерентности - использовании информации, вычисленной на предыдущем шаге симуляции движения. Алгоритм выполняет локальный поиск точек пересечения, начиная с точки, найденной на предыдущем шаге. Представив каждое тело выпуклым многогранником, можно осуществить данный поиск, переходя от вершины к вершине по ребрам и граням тела. Для корректности алгоритма Лина-Канни необходимо, чтобы многогранники были выпуклыми и непересекающимися. Алгоритм V-Clip является расширением алгоритма Лина-Канни на случай невыпуклых и пересекающихся многогранников.

Коммерческое программное обеспечение для симуляции движения

Обычно динамическая симуляция механизмов осуществляется с помощью автономных пакетов программ, относящихся к классу CAE - систем инженерного анализа. Два самых известных пакета динамической симуляции - это ADAMS (производства MSC.Software) и DADS (производства LMS; ныне DADS - часть пакета Virtual.Lab Motion). Основная функциональность таких систем позволяет:

? импортировать геометрические модели из CAD-системы;

? собрать механизм из деталей, связав их кинематическими парами;

? задать действующие силы (контакта, пружинных соединений, трения);

? вычислить массы деталей и их тензоры инерции;

? определить столкновения в процессе симуляции работы механизма;

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

Вопросы для самоконтроля

1. Какими уравнениями описывается движение системы твердых тел в поле действия сил?

2. Как моделируются контакты тел при описании динамической системы с помощью уравнений Ньютона-Эйлера?

3. Что такое уравнения Лагранжа и в чем их преимущества перед уравнениями Ньютона-Эйлера при решении задач динамики системы твердых тел?

4. Опишите общую схему методов определения столкновений. Для чего они используются в САПР?

5. Опишите алгоритмы широкой фазы при определении столкновений.

6. Опишите алгоритмы фазы сужения при определении столкновений.

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

Дополнительная литература

Освежить теоретические знания о динамике движения твердого тела поможет классический учебник [6]. Обзор известных подходов к численному моделированию динамики системы твердых тел и механизмов дан в [27]. Методы определения столкновений в системе твердых тел исчерпывающим образом охарактеризованы в [24]. Работа [4], описывающая оригинальный подход к моделированию динамики системы твердых тел, является едва ли не единственной, доступной на русском языке.



Лекция 12. Инженерный анализ методом конечных элементов

Конечно-элементный анализ

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

Конечно-элементный анализ основан на использовании математического метода конечных элементов (МКЭ).

Введение в метод конечных элементов

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

Анализ упругости тела

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

Рассмотрим трехмерный объект, который находится в равновесном состоянии под воздействием некоторой нагрузки, задаваемой:

? объемными (массовыми) силами fV, действующими на все точки объекта x х V;

? поверхностными силами (трения) fS, действующими на части поверхности объекта S;

? сосредоточенными внешними силами fi, i = 1,…, n, действующими на конкретные точки приложения xi.

Смещение точек деформируемого тела под нагрузкой (по сравнению с состоянием без нагрузки) описывается неизвестной функцией u: R3 R3, определенной для всех точек тела:

Смещения точек тела, описываемые функцией u, ведут к появлению в теле деформаций и напряжений.

Тензор деформаций

Деформацией называется изменение относительных размеров и формы тела под нагрузкой. Рассмотрим две близкие точки тела x и x + x. В состоянии под нагрузкой они сместятся в точки x + u(x) и x + x + u(x + x) соответственно. Обозначив разницу между точками под нагрузкой ?x?, получаем x = x + u(x + x) - u(x). В пределе (при x = 0) направление x (то есть dx) задается формулой dx = dx + u(x)dx = (I + u)dx, где

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

(dx1, dx2) = ((I + u)dx1, (I + u)dx2) = (dx1, (I + u)T(I + u)dx2) =

= (dx1, dx2) + (dx1, (uT + u + uTu)dx2).



Симметрическая матрица , описывающая изменение скалярного произведения между векторами, задающими направления близких точек тела, называется лагранжевым тензором деформаций. Заметим, что этот тензор, как и вектор смещений, вообще говоря, различен в каждой точке тела. В случае E* = 0 при u 0 имеем дело с бесконечно малым вращением без деформаций. При малых деформациях квадратичные члены матрицы деформаций обычно опускают, получая классический тензор деформаций:

Физический смысл коэффициентов тензора следующий. Коэффициенты ii задают линейную деформацию вдоль трех координатных осей (то есть предел отношения изменения расстояния между точками к исходному расстоянию), а коэффициенты 2ji (i j) задают угловую деформацию между координатными осями.

Тензор напряжений

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

Тогда предел существует и называется напряжением в точке P вдоль вектора n. Для определения напряжения в произвольном направлении используется тензор напряжений, который задает напряжение в произвольном направлении n как tn = n. Для большинства материалов тензор задается симметрической матрицей. Физический смысл тензора напряжений иллюстрируется на примере срезов, параллельных координатным плоскостям (рис. 31).

Рис. 31

Обобщенный закон Гука, матрицы жесткости и упругости

В предыдущей лекции мы рассматривали динамику движения недеформируемого твердого тела. Как известно, она описывается уравнениями Ньютона-Эйлера (основанными на втором законе Ньютона), которые связывают линейное и угловое ускорение тела с действующими на него силами посредством массы и тензора инерции. Аналогичную роль в задачах теории упругости играет обобщенный закон Гука. Открытый в XVII веке английским математиком Гуком (Hooke) закон растяжения для тонкого стержня имеет вид F = kx, где F - действующая на стержень сила, x - величина растяжения, а k - коэффициент упругости. Величину коэффициента упругости можно связать с физическими размерами стержня следующим образом: k = ES/L, где S - площадь поперечного сечения стержня, L - его длина, а E - модуль Юнга. С введением понятий относительного удлинения x/L и нормального напряжения в поперечном сечении = F/S закон Гука принимает вид E. Как мы уже знаем, в общем случае напряжение и деформация определяются симметрическими тензорами размера 3х3. Тем не менее, между этими тензорными сущностями действует то же линейное соотношение, что и между скалярными, называемое обобщенным законом Гука: . Тензор четвертого порядка C в данной формуле задается 81 коэффициентом (34), но так как он связывает симметрические 3х3-матрицы, каждая из которых определяется шестью скалярными величинами, то может быть представлен в виде 6х6-матрицы жесткости D, которая связывает вектор деформаций с вектором напряжений

Это линейное соотношение справедливо для малых деформаций. Обратная матрица к матрице жесткости (D-1) называется матрицей упругости. Заметим, что матрица жесткости (упругости) полностью определяется свойствами материала и не зависит от конкретных нагрузок на тело. Упругие свойства материала могут быть описаны с помощью двух параметров, измеряемых в заданных направлениях: модуля Юнга E, который определяет отношение напряжения к внутренней деформации и коэффициента Пуассона н, который характеризует отношение относительного поперечного сужения к относительному продольному удлинению . Для линейно-эластичных изотропных материалов (таких как металлы, стекло, полипропилен и полиэтилен, резина - в случае малых деформаций) модуль Юнга и коэффициент Пуассона являются константами, не зависящими от направления и точки измерения, а матрица жесткости имеет следующий вид:

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

Уравнение равновесия тела под нагрузкой

Внутренняя энергия dU бесконечно малого объема dv рассматриваемого тела под нагрузкой определяется по его напряжениям и деформациям с помощью следующей формулы:

Полная потенциальная энергия тела P = U - W равна разности внутренней потенциальной энергии тела и виртуальной работы внешних сил:

Применение МКЭ для расчета малых напряжений тела под нагрузкой

Представим рассматриваемое тело в виде совокупности конечного числа узловых точек x1, …, xn, каждой из которой сопоставим неизвестный вектор ui, определяющий значение искомой функции (смещения u) в точке xi: u(xi) = ui. Для представления значения u в других точках разобьем тело на конечные элементы с вершинами в узловых точках. Внутри каждого конечного элемента значения функции u будем интерполировать по значениям в узловых точках. Более подробно виды конечных элементов и функции интерполяции разберем ниже, а пока ограничимся простейшими тетраэдральными элементами с линейной интерполяцией. Тетраэдр Ti = (xi1, xi2, xi3, xi4) задается четырьмя узловыми точками и представляет собой трехмерный симплекс:

Коэффициенты tik называются барицентрическими координатами точки

x = ti1xi1 + ti2xi2 + ti3xi3 + ti4xi4

в тетраэдре Ti. Линейная интерполяция предполагает следующее выражение для искомой функции u в точке x, принадлежащей тетраэдру Ti:

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

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

Аналогично обстоит дело с поверхностным интегралом, который надо вычислять по граням конечных элементов, приходящихся на выделенную часть границы тела. Работа точечных сил вообще не требует интегрирования и вычисляется напрямую. В результате в качестве выражения для полной виртуальной работы W всех внешних сил получим линейное выражение от неизвестных векторов u1, …, un. Используя вектор u‚ = (u1, …, un), виртуальную работу можно задать равенством W = (f, u‚) + const.

Интегрируя приведенное выше выражение для внутренней энергии в виде суммы интегралов по конечным элементам, получаем квадратичную форму (Ku‚, u‚) + const. Согласно вариационному принципу Лагранжа вариация разности потенциальной энергии и виртуальной работы внешних сил должна быть нулевой, то есть (Ku‚, u‚) = (f, u‚) + const для всех допустимых перемещений u. Это значит, что вектор u‚ должен удовлетворять системе линейных уравнений Ku‚ = f. Таким образом, для расчета малых смещений тела под нагрузкой методом конечных элементов достаточно решить систему линейных уравнений. Решение этой системы задает смещение узловых точек сетки конечных элементов. Смещение внутренних точек сетки вычисляется затем с помощью интерполяционных полиномов. По смещениям можно вычислить тензоры деформаций и напряжений в каждой точке тела.

Другие приложения МКЭ

Так как вариационный принцип Лагранжа справедлив не только для механических систем, то метод конечных элементов находит широкое применение при моделировании самых разных физических процессов. Физическая интерпретация векторов u‚ и f дана в табл. 5. Матрицу K во всех приложениях МКЭ традиционно называют матрицей жесткости.

Таблица 5

Область приложения	Вектор состояния (u‚)	Сопряженный вектор (f)
Механика твердых тел Теплообмен Течение Электростатика Магнитостатика	Перемещение Теплопроводность Скорость Электрический потенциал Магнитный потенциал	Сила Тепловой поток Поток Плотность заряда Интенсивность магнитного поля



Типы конечных элементов

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

Рис. 32

Разбиения для МКЭ

Для автоматического построения разбиения тела на конечные элементы чаще всего применяются следующие методы:

? соединения узлов;

? топологического разбиения;

? геометрического разбиения;

? решеточные;

? отображаемых элементов.

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

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

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

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

Общая схема конечно, элементного анализа в CAE,системах

Большинство систем автоматического конструирования и инженерного анализа (CAE) позволяют проводить анализ различных характеристик изделия с помощью метода конечных элементов. На входе таких систем задается геометрическая модель изделия, подготовленная в CAD-системе. В случае интегрированных CAD/CAE-систем нет никаких сложностей с передачей модели из CAD-части в CAE, однако для автономных CAE-систем требуется перевод модели в формат, воспринимаемый данной CAE-системой. К счастью, большинство коммерческих CAE-систем понимают форматы большинства коммерческих CAD-систем, в крайнем случае, можно воспользоваться нейтральным форматов (IGES или STEP).

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

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

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



Коммерческие пакеты конечно-элементного анализа

Несмотря на то, что большинство интегрированных CAD/CAM/CAE-систем имеют встроенные модули для расчета деформаций и некоторых других характеристик изделия методом конечных элементов, по-настоящему полная функциональность для конечно-элементного анализа предлагается только в рамках автономных CAE-систем. Самыми известными из них являются:

? ANSYS (производства одноименной компании);

? NASTRAN (производитель - MSC.Software);

? ABAQUS (до недавнего времени независимая компания, ныне принадлежит Dassault Systeкmes);

? LS-DYNA (Livermore Software Technology Corporation).

На рис. 33 ниже приведен пример расчета деформаций с помощью МКЭ в системе CATIA V5.

Рис. 33

Вопросы для самоконтроля

1. Что такое конечно-элементный анализ? На каком математическом аппарате он основан? Каковы области его применения?

2. Опишите постановку задачи расчета деформации тела под нагрузкой.

3. Что такое тензоры деформаций и напряжений? Охарактеризуйте их физически и математически.

4. Опишите обобщенный закон Гука. Какие свойства материала определяются модулем Юнга и коэффициентом Пуассона?

5. Опишите и объясните уравнение деформации тела под нагрузкой.

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

7. Что такое обобщенная матрица жесткости? Опишите применение метода конечных элементов для различных классов физических задач.

8. Опишите основные типы конечных элементов.

9. Какие существуют способы построения сеток для метода конечных элементов?

10. Опишите общую схему конечно-элементного анализа в CAE-системах и приведите примеры пакетов программ.

Дополнительная литература

Краткое введение в конечно-элементный анализ можно найти в восьмой главе книги [8], а также в третьей части [5]. Более подробное изложение содержится в монографии [6] и в лекциях [9].



Лекция 13. Автоматизация производства

Архитектура станков с ЧПУ

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

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

Станок с ЧПУ (numerical control, NC) отличается от обычного наличием блока управления станком (machine control unit, MCU), функционально состоящего из модуля обработки данных (data processing unit, DCU) и замкнутой системы управления (control loop unit, CLU). Модуль MCU считывает данные с входного источника данных (сначала это были перфокарты и перфоленты, потом их заменили магнитные носители), а DCU преобразует их в сигналы управления станком. С появлением микросхем и микропроцессоров эта архитектура эволюционировала в привычную современную схему управления любым внешним устройством с помощью интегрированного в устройство контроллера и внешнего компьютера. С введением компьютера в схему управления станком последний стал называться станком с компьютеризированным числовым управлением (computer numerical control, CNC). В современных производственных цехах все компьютеры, контролирующие станки с ЧПУ, соединены в сеть под командой центрального компьютера, с которого и происходит непосредственное управление всем цехом, включая загрузку данных на конкретный станок. Подобная схема называется распределенным чис- ловым управлением (distributed numerical control, DNC).

Принципы программирования для станков с ЧПУ

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

? N - идентификатор блока;

? G - предварительные команды (позиционирование, остановка, ускорение и пр.; все доступные команды закодированы двузначными числами);

? X, Y, Z, I, J, K - координаты инструмента в системе координат станка (X, Y, Z - координаты конечной точки, I, J, K - координаты центра окружности, задаются только при движении по дуге);

? F - скорость подачи (в дюймах в минуту);

? S - скорость шпинделя (в оборотах в минуту);

? T - номер инструмента (резца, фрезы, сверла - зависит от станка);

? M - вспомогательные команды (включить охлаждение и пр.;

все доступные команды закодированы двузначными числами).

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

N009 G02 X-28.28 Y0.0 I14.14 J5.0

Семантика команды такова: из текущей точки резец должен описать дугу окружности с центром в точке (14.14, 5.0) в направлении движения часовой стрелки и остановиться в точке (-28.28, 0.0). Если после этого надо с повышенной скоростью (не касаясь обрабатываемой детали) перевести резец в точку (-75.0, 0.0, 40.0), после чего выключить шпиндель и охлаждение, то соответствующая команда будет выглядеть так:

N010 G01 X-75.0 Y0.0 Z40.0 F950 M30

Система координат целиком определяется типом станка. Существует общее соглашение, что ось z совпадает с осью вращения шпинделя (на котором закреплена деталь или инструмент), а ось x всегда горизонтальна. Причем движение в положительном направлении по обеим осям удаляет инструмент от заготовки.