CC BY
32
Серия 3. Естественные науки.
Алгоритмические концепции метода твердых тел
к.т.н. доц. Русанов П.Г. Университет машиностроения 8 (495) 223-05-23 доб. 13-18, [email protected]
Аннотация. Метод твердых тел входит в группу адаптивных вычислительных методов механики сплошной среды. Он применим для анализа движения твердых тел, жидкостей и газов. Конструктивными приемами метод обеспечивает априорное разделение обобщенных координат для конечных объемов объекта исследования на быстрые и медленные переменные. Общее число медленных переменных не превышает 6 N, где N - число конечных объемов. Указаны способы формирования математической модели состояния объекта относительно медленных переменных без участия быстрых переменных.
Ключевые слова: твердое тело, жидкость, газ, дискретный элемент, собственный репер, глобальные и локальные переменные, модели инерционных и физических сил
Введение
Исторически теоретические образы вещества «первого» приближения в механике в виде материальной точки и абсолютно твердого тела (АТТ) послужили отправными катализаторами научного и технического прогресса человеческой деятельности. Затем к ним примкнули модели деформируемого твердого тела, жидкости и газа, для поведения анализа которых потребовались специальные методы. На рынке компьютерных услуг лидерами среди них в настоящее время стали метод конечных элементов (МКЭ) и метод конечных разностей (МКР), каждый из которых «с радостью» переложил проблемы аналитической разрешимости уравнений состояния на «плечи» «так удачно подвернувшейся» вычислительной техники, правда, ценой потери наглядности исходного физического образа исследуемого объекта. В стратегии указанных методов заложены различные способы математической аппроксимации искомых функций их значениями на дискретном множестве точек пространства и времени, называемом сеткой математической дискретизации. Поэтому ключевыми элементами технологий сеточных методов стали алгоритмы численной реализации математических преобразований. Для получения более точных результатов эти методы рекомендуют применять более мелкую сетку дискретизации, которая, к сожалению, слишком «жадно поглощает» ресурсы памяти компьютерной техники.
В данной статье рассмотрен иной концептуальный подход к решению проблемы численного моделирования поведения механических объектов, названный методом твердых тел (МТТ). В отличие от сеточных методов, МТТ приходит к единой методологии формирования математической модели (ММ) для всех типов веществ с укрупненной сеткой дискретизации массы объекта исследования, физически близкой к образам АТТ. Сочетание простых алгоритмов формирования ММ с высокой точностью результатов позволяет методу МТТ конкурировать с сеточными методами в задачах с преобладанием механических форм движения вещества.
Метод опирается на ряд конструктивных приемов и аналитических разработок:
• условное разбиение объекта исследования на дискретные элементы (ДЭ);
• конструктивное введение обобщенных координат для системы ДЭ, априорно состоящих из медленных и быстрых переменных;
• предварительная аналитическая оценка инерционных и физических сил элементов с помощью медленных переменных, без участия быстрых переменных.
Среди достоинств метода МТТ отметим следующие:
• метод применим для анализа механического движения всех видов вещества;
• объёмный ДЭ имеет не более 6-ти степеней свободы при минимуме, равном 12, у объем-
ного конечного элемента в методе МКЭ;
• МТТ располагает уникально простыми аналитическими оценками инерционных и физических сил дискретных элементов;
• МТТ сводит задачу исследования к анализу поведения лишь одних медленных переменных;
• в зависимости от постановки задачи ММ имеет вид системы алгебраических уравнений или нежесткой системы обыкновенных дифференциальных уравнений.
Условные обозначения и исходные ограничения
Для упрощения изложения основ метода вынужденно ограничим постановку проблемы. В качестве объекта исследования рассматриваем одиночную однородную среду, которая на протяжении всего времени наблюдения за ней допускает условное разграничение на N дискретных элементов (ДЭ), сохраняющих свои массы. В случае твердого тела его внешнюю и внутреннюю естественные поверхности также относим к сетке дискретизации. Для жидкости и газа номинальную феноменологическую форму поверхности ДЭ считаем сферической.
Во избежание неоднозначности в вопросе принадлежности материальных точек к тем или иным ДЭ будем считать, что у разных элементов общими могут быть только безмассовые точки сетки дискретизации. Пока два или более элементов имеют хотя бы одну общую точку, они считаются соседними.
Для твердого тела общие участки поверхности сетки дискретизации двух соседних ДЭ назовем их условными границами, кратко обозначаемыми далее как «У-границы». А участок сетки дискретизации ДЭ, принадлежащей естественной поверхности твердого тела, будем называть его естественной границей и обозначать кратко как «Е-граница».
Непринципиально сузим постановку динамической задачи тем, что первоочередной целью ее решения будем считать расчет кинематических характеристик пространственного движения объекта в заданных физических условиях, с известными начальными условиями на заданном интервале времени. Очевидно, что на основании такого расчета далее может быть получена информация и о силовых взаимодействиях. Поэтому сосредоточим внимание только на способах формирования ММ, замкнутой относительно кинематических переменных и не содержащей характеристик пассивных сил. В связи с этим наложенные на среду связи считаем идеальными. Для сужения темы пока пренебрежем тепловым движением масс, а также влиянием электромагнитных полей и излучений.
Неизвестное поведение сплошной среды будем изучать в глобальной системе координат (ГСК), которую временно лишим самостоятельного движения в инерциальной системе координат, то есть считаем её также инерциальной.
Обобщенные координаты для ДЭ. Первично каждый ДЭ рассматриваем как систему из п индивидуально движущихся точечных масс, обладающих Ъп степенями свободы относительно ГСК. Воспользуемся классическими обозначениями, опустив порядковый номер ДЭ
среди других: ^ = ~ °бщая масса ДЭ, тмасса к-той точки, точка С - его центр масс
1
(ЦМ), / - матрица инерции масс ДЭ в центральных осях. Для краткости изложения понятие «главные центральные оси инерции массы» ДЭ заменим термином «собственный репер» (CP). Геометрический образ CP будет выполнять важную, конструктивную роль в обосновании и выводах метода.
Чтобы «сблизить» расчетные образы одного ДЭ и некоторого АТТ той же массы, специально «внедрим» в число обобщенных координат (ОК) материальных точек ДЭ 6 новых переменных, традиционно используемых в расчетах динамики АТТ. Для этого тривиальный список из Ъп индивидуальных координат каждой материальной точки в ГСК заменим другим, состоящим из 2-х групп кинематических параметров - глобальных и локальных перемен-
ных (ЛП). К глобальным переменным (ГП) отнесем 6 параметров, задающих положение собственного репера в ГСК в текущий момент времени. Причем для определенности далее будем считать, что в качестве первых трех ГП, выбраны Хс, Ус, три декартовые координаты ЦМ, а остальные три ГП - это три угла пространственной ориентации собственного репера ДЭ в осях ГСК, которые условно обозначим как а, Р,ф .
Остальными обобщенными координатами, относимыми к группе локальных переменных, будут служить Зп - 6 независимых декартовых координат материальных точек ДЭ в осях СР с учетом составного способа задания их положения по отношению к ГСК,
гок = гс + А\, (к= 1, п),
где: гс = (Хс,УС,Zc), гок = (Хк,Ук,2к) - векторы положения ЦМ элемента и к- той
материальной точки ДЭ в ГСК; Тк = (хк,ук, 1к) - вектор положения точки в СР; А -
матрица преобразования поворота от осей ГСК к осям СР, задаваемая углами а, Р,ф. Не имея возможности в явном виде указать на 6 зависимых из общего числа 3п декартовых координат точек в СР, воспользуемся чисто математическим способом их выявления на основании уравнений связи, отражающих инвариантные геометрические свойства расположения точечных масс ДЭ в осях СР для любого момента времени существования ДЭ:
• равенство нулю трех статических моментов точечных масс ДЭ относительно координатных плоскостей СР
п п п
хк=Ук=Итк 2к=°> (!)
1 1 1
• равенство нулю трех центробежных моментов инерции точечных масс в осях СР
п п п
Т,тк хкУк = Т,тк хкЧ =Е тк УкЧ = 0' (2)
1 1 1
Предпринятое конструктивное расслоение ОК позволит далее обосновать переход от формальной ММ в виде набора из 3пЫ дифференциальных уравнений (ДУ) для точечных масс всего объекта к более компактной подгруппе ДУ, замкнутой относительно 6Ы глобальных переменных системы собственных реперов ДЭ. Новую форму ММ для объекта исследования с участием одних ГП будем называть его генерализованной моделью (ГМ).
По терминологии метода разделения быстрых и медленных движений МТТ отводит глобальным переменным роль медленных переменных. Они будут служить основными информационными неизвестными задачи исследования. Метод рассматривает ЛП как быстрые переменные, которые, однажды возникнув, по ряду физических причин достаточно быстро сходятся к своим постоянным значениям за доли самого короткого из периодов собственных колебаний всей системы СР.
В связи преднамеренным отказом метода от изучения поведения огромного количества ЛП, ради значительного упрощения всей проблемы исследования, сама идеология метода будет выполнять роль своеобразного низкочастотного фильтра, «не пропускающего» в ГМ высокочастотные движения («флюктуации») точечных масс. При этом метод ставит задачу пренебрежения лишь флюктуациями ЛП, а не самими величинами ЛП.
Блоки математической модели. Проследим за «нежелательным» проникновением флюктуаций ЛП в генерализованную модель. От каждого ДЭ включим в ГМ по 6 независимых уравнений динамики, образующих 2 блока:
• блок 1 - это три ДУ на основании теоремы о движении ЦМ:
Ма = К, _ (3)
• блок 2 - это три ДУ на основании теоремы об изменении К - главного вектора кинетического момента вокруг ЦМ от относительных количеств движения точечных масс в осях
Серия 3. Естественные науки.
Кенига (за оси Кенига принимаем постоянную во времени копию направлений осей ГСК, но с началом в подвижной точке С):
К = Ь, _ _ (4)
где: а - вектор ускорения ЦМ ДЭ относительно ГСК, Я, Ь - главные векторы внешних для ДЭ физических сил и моментов этих же сил вокруг ЦМ. Поочередно рассмотрим проявления «флюктуаций» ЛП, сначала в инерционных членах
ГМ, а затем и в силовых характеристиках Ь .
Модель инерционных сил ДЭ. Не без некоторой доли морального удовлетворения отмечаем, что теорема (3) тоже частично «придерживается» стратегии изоляции ЛП в инерционных членах, поскольку вектор а определяется вторыми производными по времени только от первых трех назначенных нами ГП и поэтому принципиально не зависит от ЛП. Менее благоприятна ситуация в (4), где вектор К зависит от ОК обоих типов:
п
К = хтк(шхгк + щ)] = Ке + Кг, (5)
1
где: гк, йк = Тк, СО, К, Ке, Кг - векторы, представленные своими проекциями на оси СР; Гк, ик векторы положения и скорости массы тк относительно СР;
Ю - вектор угловой скорости сферического движения СР относительно осей Кенига, или ГСК. Проекции вектора Ю на оси СР определяются кинематическими формулами Эйлера - Пуассона для комбинации трех углов, выбранных в качестве глобальных переменных;
п _
Ке =^ \гк х тк (шх 7к)] = 1ш - часть вектора К от переносных скоростей точечных
1
масс ДЭ вследствие сферического движения СР относительно осей Кенига; / - диагональная матрица моментов инерции масс ДЭ в осях СР;
п _
Кг = ^ (гк х ткйк) - часть К от ик - относительных скоростей точек ДЭ в осях СР.
1
Считая, что векторы, входящие в (5), заданы проекциями на подвижные оси СР, вычисляем абсолютную производную по времени от вектора К по формуле:
К = К + шх К, (6)
где: К - вектор локальной скорости изменения вектора К по отношению к осям СР. На основании (5) и (6) получаем:
К = Ке +Кг; Ке = Ке + шх Ке; Кг = Кг + шх Кг,
Ке = 7Ш+ /Ш+ шх 7ш, (7)
. п п
Кг = хтк ) + (гк хтк йк) (8)
1 1
Заметим, что в частных случаях последнее слагаемое в (7) может тождественно равняться 0, если вектор Ю все время направлен по главной оси инерции ДЭ. Так будет, например, у элемента со сферическим эллипсоидом инерции.
Из формул (7), (8) следует, что сами ЛП и их производные входят во все инерционные члены уравнения (4). И, несмотря на предполагаемую ограниченность диапазонов изменения величин самих переменных х^), Ук(1), ¿¿(0 и «обнадеживающие» следствия из (1) и (2) в виде:
га га
1 1
га га га
Т,тк (хкУк +хкУк) = Т,тк (хкЧ + хк ) = Xтк (УкЧ + Ук ) = О, (10)
1 1 1 тем не менее, в этой ситуации отсутствуют чисто математические основания для освобождения выражений (7) и (8) от производных ЛП.
В связи с этим делаем два вывода.
Вывод 1. Из условий (1), (2), (9) и (10), которым подчинено поведение ЛП, сами по себе не являются достаточными для пренебрежения их величинами в инерционных членах (7), (8) уравнения (4). Для чисто формального исключения «флюктуаций» ЛП из левой части (4), достаточно одновременного выполнения двух, в общем случае независимых, условий: а) элементы матрицы /(- постоянны и б) вектор Кг (0 - тождественно равен 0, которые, очевидно, можно заменить одним, но еще более строгим условием:
ь (0 = А (0 (0=0(*=м )• (11)
Условие (11) соответствует пренебрежению кинетической энергией движения точечных масс ДЭ относительно СР или, что то же самое, принятию допущения о постоянстве положений точечных масс в объеме каждого ДЭ. В таком случае, помимо изначально постоянной массы ДЭ, в осях СР будет постоянной во времени и матрица /, и расчет 1С можно
выполнять по формуле, справедливой для АТТ: К = ЛаН-шх 1ш.
Корректность принятия допущения (11) в отношении твердых тел можно аргументировать с помощью опытных данных. Например, при частых ударах кончиками пальцев по корпусу гитары четко различимы на слух отдельные звуки после каждого удара. То есть ударные возмущения корпуса затухают быстрее, чем наносится очередной удар. Тем самым условие (11) можно считать вполне приемлемым на конечных интервалах времени, соответствующих времени релаксации для низших гармоник.
Кроме того опыт свидетельствует о том, что время затухания собственных колебаний в теле резко сокращается с увеличением их частоты. Похоже, что это явление можно объяснить наличием «резонансной откачки» энергии у возбужденных участков тела и передачи ее на соседние области тела в зоне близких собственных частот. Неслучайно исследователи колебаний механических систем, зная о скоротечности затухания высокочастотных форм колебаний, проявляют свой практический интерес лишь в отношении низших собственных частот колебаний.
В теории колебаний установлено, что частоты собственных колебаний упругого стержня обратно пропорциональны квадрату его длины. Поэтому сама процедура «дробления» объекта исследования на ДЭ объективно приводит к смещению спектра собственных частот колебаний ДЭ в область повышенных значений. Приведенные рассуждения и факты отражают асимптотический характер снижения погрешности получаемых результатов на основании метода МТТ. Очевидно, что при измельчении сетки дискретизации, направленном на увеличение N - числа ДЭ при одновременном уменьшении их объемов, менее «подозрительным» станет допущение (11) и само решение задачи будет более точным, так как в него войдут еще более высокие гармоники собственных колебаний.
Условие (11) не эквивалентно известному допущению о «неизменности начальных размеров», применяемому в отношении достаточно жестких конструкций. С этим замечанием легко согласиться, если вспомнить неразрушающий механический прием, позволяющий отличить куриное яйцо, сваренное «вкрутую», от неварёных яиц. Т.е. сохранение формы объекта еще не гарантирует однозначности его поведения в опыте.
Вывод 2. Допущение (11) «освобождает» от «флюктуаций» ЛП инерционные члены в (3) и (4), и поэтому ГМ не будет иметь сингулярных возмущений. Оно же позволяет нам качественно изменить отношение к изучаемому объекту и далее его рассматривать уже не как систему материальных точек., а как систему ДЭ с постоянными массами и постоянными центральными осевыми моментами инерции массы, преемственно отражающую исходные инерционные свойства объекта исследования.
С научной точки зрения, в принятом допущении (11) нет как такового элемента новизны, кроме необычности его применения в отношении отдельного ДЭ. Науке известно большое число примеров, где весь объект исследования (целиком), невзирая на наличие у него «податливых» мест, считался абсолютно твердым. Поэтому, следуя историческому опыту моделирования механических объектов, вполне естественно применять и более гибкую разновидность расчетной схемы, в которой не только все тело целиком, но и его части также претендуют на обладание инерционными свойствами абсолютно твердых тел.
Очевидно, что и модель одиночного АТТ, и модель континуума становятся частными случаями расчетных схем метода МТТ, поскольку при N= 1 ДЭ переходит в разряд одиночного АТТ, а при неограниченном росте числа ДЭ и стремлении объема каждого ДЭ к нулю МТТ приходит к дискретно-точечной расчетной схеме среды.
Замечание. Методы расчета положения ЦМ, величины массы, направлений осей СР и величин моментов инерции ДЭ по заданной форме и распределению плотности материала не имеют существенных математических проблем и поэтому здесь не рассматриваются.
Внешние физические силы ДЭ. Теперь проследим за участием ГП и ЛП в выражениях
для R,L- главных векторов внешних сил и моментов внешних сил вокруг ЦМ ДЭ системы (3),(4). Одновременно с этим будем искать возможные варианты оценок этих векторов лишь при помощи ГП и без участия ЛП. Вначале сделаем ряд важных замечаний.
В механике силовые взаимодействия тел подразделяют на распределенные по поверхности или по объему, называя их соответственно силами близкодействия и дальнодействия. Поэтому употребление термина «сила» далее будет справедливым лишь в отношении понятий «главный вектор сил» и «главный момент сил», по которым оценивается уровень силового воздействия на участок поверхности или на элемент объёма.
Точность моделей силовых взаимодействий в целом невысока и находится пока на уровне погрешностей в задании физических величин, таких как коэффициент трения скольжения, механические константы упругости материала, а также в оценках сил сопротивления воздушной среды и потерь на внутреннее трение. Отсюда можно сделать вывод, что «не стоит ломать копий» в описании частных силовых моделей, если в них низка точность физических констант или одновременно с ними применяются другие, заведомо «несовершенные» силовые модели.
Последующее изложение в большей степени будет посвящено проблеме силовых моделей для деформируемого твердого тела, поскольку для газа или жидкости выбор сферической формы ДЭ практически снимает остроту изоляции ЛП в силовой модели при идеальных связях. В самом деле, при пренебрежении силами трения в контактах L - главный вектор моментов в (4) от нормальных сил для ДЭ сферической формы тождественно равен 0, что позволяет далее исследовать лишь движение ЦМ дискретного элемента .
Модель внешних физических сил дальнодействия. В рамках принятой постановки задачи объемные силы дальнодействия представлены силами гравитации. В пределах малого объема тела распределение сил гравитации от внешнего удаленного источника большой массы принято считать однородным. Поэтому главный момент таких сил вокруг ЦМ ДЭ равен нулю, а главный вектор зависит от положения ЦМ в поле гравитации, т.е. в конечном счете от ГП и без какого бы то ни было участия ЛП. Данную силовую модель универсально применяют для всех форм ДЭ и всех типов веществ: газа, жидкости, деформируемого твердого тела. Вблизи поверхности планеты для жестких тел малой массы локальным гравитационным
взаимодействием между их точечными массами принято пренебрегать.
Модели внешних физических сил близкодействия. В отличие от сил дальнодействия силовые модели близкодействия у разных веществ (газ, жидкость, твердое тело) не обладают подобным единообразием вследствие принципиальных различий в поведении структурированных и неструктурированных материалов.
Силовые модели для жидкости и газа. Феноменологические силовые модели взаимодействия дискретных элементов ДЭ из неструктурированных веществ (жидкость и газ), представляемых в МТТ в виде системы «скользких» сферических масс, могут наследовать формулы теории упругого удара, теории контактных деформаций, теории поверхностных пленок, допускающих их запись лишь через ГП центров масс ДЭ.
Естественно, что МТТ гарантирует повышение точности расчетов с увеличением числа сферических ДЭ в области с повышенным градиентом скоростей ЦМ элементов.
Силовые модели для твердых тел
У одного ДЭ из структурированного материала (случай деформируемого твердого тела) возможны три варианта реализации силового близкодействия:
• контакт ДЭ с жидкостью или газом по Е-границе;
• контакт ДЭ одного тела с другим твердым телом по Е-границе;
• взаимодействие двух ДЭ одного тела по У-границам.
Для силовых взаимодействий первого типа в большинстве случаев наука не располагает исчерпывающе точной силовой моделью. Поэтому исследователи оценивают силовое взаимодействие на основании приближенных подходов с применением феноменологических моделей, в которые входят наиболее «доступные» физические величины: площадь свободной поверхности и закон распределения векторов давления по этой площади. При учете скоростей движения частиц жидкости или газа «относительно» твердого тела его вынужденно принимают за АТТ. Причем в этих моделях более «неопределенными», в первую очередь, считаются векторы скоростей малых частиц жидкости или газа, а не вектор скорости микроучастка свободной поверхности ДЭ. Поэтому в условиях «неточных» исходных данных или пренебрежений другими источниками сил пока остаются вполне приемлемыми упрощённые силовые модели, не учитывающие вибрации отдельных частей твердого тела.
Подобная ситуация имеет место и в отношении силовых моделей контактного взаимодействия твердых тел. Пока наука находится в стадии поиска таких моделей, которые способны совместно учитывать ударный характер взаимодействия тел с энергетическими потерями на трение и износ в контакте. Как правило, применяемые модели строятся на соотношениях теории контактных деформаций для статических условий, в которых величина сближения тел определяется по точкам контр-тел, удаленных от зоны контакта, а скорости точек на поверхности тел в зоне контакта рассчитываются по формулам кинематики АТТ. Именно поэтому для оценки контактных взаимодействий ДЭ с контр-телом по таким моделям «пока» будет достаточно учитывать в них лишь поведение его глобальных переменных.
В целях унификации силовых ГУ на Е-границах ДЭ действующие на них распределенные силы близкодействия 1-го и 2-го типов представим псевдоэквивалентной системой сил, приведенной к ЦМ этого ДЭ и оцениваемой по формулам теоремы Пуансо для АТТ. В связи с этим далее будем считать, что Е-границы всех ДЭ свободны от силовых воздействий, т.е. ГУ в точках Е-границ отвечают нулевые значения напряжений.
Модель упругого взаимодействия двух ДЭ. Перейдем к третьему случаю силового близкодействия по У-границе между двумя соседними элементами дискретизации, один из которых назовем опорным, а другой - смежным. Покажем, что силовое воздействие на опорный ДЭ со стороны смежного ДЭ можно оценить в аналитической форме по величинам глобальных переменных, задающих положения СР0, СР - их собственных реперов в ГСК.
В целях упрощения изложения методики расчета полагаем, что:
• материал тела является однородным, изотропным и линейно упругим; 130 Известия МГТУ «МАМИ» № 3(17), 2013, т. 1
• в любой точке ДЭ модули элементов матрицы деформаций малы в сравнении с 1;
• физическая модель деформирования материала не относится к числу моделей моментной теории упругости;
• в недеформированном состоянии тела все его ДЭ имеют форму одинаковых кубиков, а направления их осей СР совпадают с направлениями соответствующих ребер ДЭ.
Следствия, вытекающие из принятых допущений:
• о линейном характере зависимости приведенных сил от малых независимых смещений СР смежного ДЭ относительно СРо опорного ДЭ;
• о допустимости оценки распределенных сил на У-границах без учета их деформаций, аналогично гипотезе Бернулли в теории изгиба стержней;
• о допустимости применения принципа суперпозиции для расчета сил.
С учетом принятых допущений индивидуальное силовое воздействие смежного ДЭ на опорный ДЭ будет определяться линейной зависимостью
F=C 6. (12)
Здесь Р=(ДХ, Яу, Ьх, Ьу, Ьг)т - математический вектор приведенных сил воздействия на опорный ДЭ, обусловленных малыми перемещениями СР смежного ДЭ относительно СР0; 5 = (и, V, w, а, (3, ф)т -математический вектор малых параметров, задающих текущее положение СР смежного ДЭ относительно СР0; С - постоянная матрица жесткости; (Ит Р = (Ит 8 = 6, (Ит С = 6x6.
Переменные и, V, ^ - имеют смысл малых линейных смещений ЦМ смежного ДЭМ, соответственно, вдоль осей х, у, г СР0, а а, (3 , ф - являются малыми углами поворотов от осей СР0 к осям СР. Причем в недеформированном состоянии материала 8 = (0,0,0,0,0,0) т. Выразим элементы вектора 8 через ГП собственных реперов СР0, СР.
Формулы связи кинематических параметров
Чтобы выразить а, (3, ф череза0, ро, фо и а, Р, ф, воспользуемся тем, что направления осей СР в ГСК можно задать двумя способами: 1) непосредственно, с помощью матрицы А(а, Р, ф) и 2) составным способом с применением матриц В(а, (3 , ф), Л0(а0, Ро, Фо), гдеВ(а , |3, ф) - матрица преобразования поворота от осей СР0 опорного ДЭ к осям СР смежного ДЭ. Тогда:
А(а, р, ф)=В(а, р, ф)»А0(а0, р0, Фо). (13)
С учетом свойств матриц направляющих косинусов перепишем (13) в виде:
В = А ■ Лт- (14)
Зависимости матриц А0, А, В от соответствующих углов Крылова структурно одинаковы. Для матрицы А она имеет вид
ап а12 а13 ' сР Сф от 5ф + sa ^Р сф sa 5ф - от sp сфЛ
А = ахх й22 а23 = -ф 5ф са сф - sa ф 5ф ОТ Сф + ОТ sp 5ф ' (15)
\ а\\ аЪ2 аъъ V -sa сР ОТ сР у
здесь для краткости записи тригонометрических функций sin и cos от углов а, Р, ф вместо самих функций указаны лишь первые буквы в их названиях.
Для малых осевых деформаций материала (порядка в « 0,001) абсолютные величины и/2а, v/2a, w/2a, а, (3 , ф - имеют тот же порядок малости. При уровне точности 10% для В не важен порядок поворотов на углы а, |3 , ф и допустимы приближенные оценки величин:
sin а, = (1, sin (3 = (3, sin ф = ф, cos (l=cos|3=cos ф = 1.
Тогда матрица В примет вид:
В =
1
-Ф
Ф 1
-Р
(16)
(17)
а
.Р -а
и с учетом (16)и (14) получим, что
а = ¿23, (3 = ¿31, Ф = ¿12,
где: ¿23, ¿31, ¿12 - элементы матрицы, равной произведению двух матриц А • А^.
Теперь выразим и, V, - линейные перемещения через ГП на основании векторного
уравнения (18), отражающего взаимосвязь расположений ЦМ элементов (рисунок 1).
---------
Рисунок 1. Условная схема взаимных расположений систем координат ГСК, СР0 и СР
В осях ГСК векторы положений центров масс двух соседних ДЭ связаны соотношени-
Г = Го + АТ(Р + (18)
Здесь г о = (Х0,У0,Г = (X, 7, - векторы положения ЦМ в ГСК; р = (р^ ,ру, рг)
- вектор положения ЦМ смежного ДЭ относительно СР0 в недеформированном состоянии материала, заданный проекциями на оси СР0; £ = (и, V, М?) - вектор деформационного перемещения ЦМ смежного ДЭ относительно СР0, заданный проекциями на оси СР0. Разрешим (18) относительно вектора £
£ = А0(г - Го)-р. (19)
На основании (19) получаем зависимости для и, V, ч*- локальных параметров положения СР смежного ДЭ в СР0 от ГПХ0, 70, 2о, ао, Ро, Фо, X, 7, 2, а, Р, ф.
и= (Х„ -Х0 ) ап + (Г„ -Г0 ) а21 + (2п -Ъ{)) аЪ\- рх,
у= (Х„ -Хо ) Д12 + (У„ -Уо ) а22 + {2„ -Ъ{)) а32 - ру, (20)
^ = (Х„ -Х0 ) «13 + (Г„ -Г0 ) а23 + (2п -Ъо ) аъъ- рг. Здесь щ— элементы матрицы Ао, а /,7 - индексы ее строк и столбцов. При этом формулы (17), (20) указывают на нелинейный характер участия ГП в векторе 5 силовой модели (12).
На завершающей стадии формирования силовой модели (12) рассмотрим алгоритм аналитического расчета постоянных коэффициентов матрицы жесткости С. Расчет значений этих коэффициентов производим методами тарировки на основании принципа суперпозиции. Задаваясь поочередно ненулевым значением каждой из 6-ти величин и, V, а, (3, ф и считая остальные равными 0, производим расчет Г=(КХ, Ьх, Ьу, Ьг)Т -вектора приведен-
ных сил методами теории упругости. Тогда, согласно (12), Су = Рг■/ 5/, где у - номер ненуле-
вой компоненты вектора 8. Алгоритм расчета описан ниже.
Данный подход к расчету элементов матрицы жесткости идейно близок к методике расчета приведенных сил для узлов конечных элементов в методе МКЭ. МТТ представляет его в виде ряда достаточно прозрачных, процедурных этапов:
1. Введение И-области - области ненулевых деформационных перемещений между собственными реперами соседних ДЭ.
2. Расчет габаритов участков на У-границах опорного ДЭ, оказавшихся в И-области.
3. Организация цикла из 6-ти независимых ненулевых вариантов значений и, V, ж, а, |3 , ф -относительных смещений СР смежного ДЭ по отношению к СР0. На каждом шаге этого цикла выполняются следующие процедурные действия:
3.1. Генерация в аналитическом виде законов деформационных перемещений материала внутри И-области, удовлетворяющих ГУ в фиксированных точках на границе области, а также в ЦМ смежного ДЭ.
3.2. Расчет по формулам теории упругости (кинематические соотношения Коши и обобщенный закон Гука) законов распределения напряжений в И-области.
3.3. Внутренний цикл для каждой У-границы опорного ДЭ, оказавшейся в И-области:
3.3.1. Аналитический расчет для опорного ДЭ главного вектора и главного момента сил, распределенных на его У-границе с приведением к одной из ее точек - центру грани.
3.3.2. Перерасчет главного вектора и главного момента сил к ЦМ опорного ДЭ.
3.3.3. Суммирование результатов для главного вектора и главного момента сил, приведенных к ЦМ опорного ДЭ с аналогичными результатами для других У-границ
ДЭ.
4. Расчеты коэффициентов в столбце матрицы жесткости, номер которого равен номеру ненулевой компоненты в векторе 8, по результатам, полученным в пункте 3.3.3.
Пояснения к алгоритму
1. Считается, что вне И-области деформации материала отсутствуют. Причем И- область может быть двух видов и Б2. Это связано с типами соседних ДЭ. ДЭ будем считать внутренним, если у него нет Е-границ, в иных случаях, ДЭ считаем внешним. Необходимость введения разграничения ДЭ на 2 типа связана с различием требований МТТ к формулировке граничных условий для законов распределения деформационных перемещений и напряжений. Для пары соседних ДЭ, из которых хотя бы один является внутренним, расчетная область имеет тип а в иных сочетаниях соседних ДЭ - область Б2.
Ограниченная ^-область представляет собой выпуклый объём, поверхность которого проходит через ЦМ опорного ДЭ и ЦМ всех ДЭ, являющихся одновременно соседними как для опорного, так и для смежного элементов.
Ограниченная .Ог-область, как и ^-область, тоже представляет собой выпуклый объём, одна часть поверхности которого проходит через ЦМ опорного ДЭ и через ЦМ всех ДЭ, являющихся одновременно соседними как для опорного, так и для смежного элементов, а другая часть совмещена с Е-границами опорного и соседних с ним ДЭ.
Очевидно, что для упрощения расчета границ участка интегрирования распределенных сил области .Ог удобно очерчивать поверхностями первого или второго порядков.
2. В ^-области поле деформационных перемещений должно удовлетворять только кинематическим ГУ в точках ЦМ, расположенных на границе ^-области соседних ДЭ, независимо от их типов.
В П2-области поле деформационных перемещений должно удовлетворять кинематическим ГУ в точках ЦМ, расположенных на границе П2-области, и нулевым значениям нормальных и касательных напряжений в центральных точках Е границ.
3. При наличии элементов симметрии или косо-симметрии в распределении напряжений на У-границе можно ограничиться расчетом лишь тех приведенных силовых факторов,
которые отвечают данному виду симметрии, а остальные принять равными 0.
4. Для регулярной сетки дискретизации в виде одинаковых кубиков в недеформированном состоянии материала возможны лишь три топологически различных варианта расположения смежного ДЭ по отношению к опорному ДЭ. Этот факт позволяет математически трансформировать однажды полученные матрицы жесткости для трех базовых случаев соседства ДЭ на все прочие случаи расположений СР в пространстве в СР0.
5. Если в недеформированном состоянии тела из однородного материала ЦМ всех ДЭ расположены на одной регулярной линии (одномерный случай) или на одной регулярной поверхности (двумерный случай), то формуле (12) будут отвечать формулы сопротивления материалов и строительной механики для стержней и пластин.
6. Один из способов существенного снижения затрат при расчете силового взаимодействия основан на общепринятом пренебрежении некоторыми видами деформаций, например для стержня и пластин, растяжением продольной оси стержня или деформациями сдвига в поперечных сечениях, что равносильно введению искусственных узлов кинематических соединений между ДЭ. В подобных случаях может полностью отпадать необходимость расчета приведенных сил и моментов, поскольку их величины определены в явном аналитическом виде в инженерных справочниках расчета на прочность.
Аналитические расчеты матриц жесткости С для базовых вариантов взаимного расположения пары соседних ДЭ несложно провести с помощью компьютерных пакетов символьных вычислений, например, таких как пакет Мар1е.
Отметим, что для двух соседних ДЭ можно оценивать силовое взаимодействии не во всей И-области, а лишь в пределах одной общей У-границы на основании более «грубой» силовой модели взаимодействия ДЭ, использующей «стержневую аналогию». То есть, приведенные силы воздействия смежного ДЭ на опорный ДЭ на У-границе можно рассчитывать на основании формул сопротивления материалов и строительной механики для различных вариантов взаимного расположения опор стержня. Здесь существенно то, что помимо сил дальнодействия суммарное воздействие на опорный ДЭ будет определяться по всем смежным для него ДЭ.
Проведение компьютерного анализа динамики стержней и нитей на основе метода МТТ доступно студентам ВУЗа начиная со 2 - 3-го курсов обучения, после усвоения классических курсов теоретической механики и сопротивления материалов.
Технология построения математических моделей по методу твердых тел прошла практическую апробацию на ряде численно решенных задач и подтвердила свою эффективность в отношении таких проблем, как анализ статических форм, устойчивости, форм и частот собственных колебаний, а также динамических процессов для различных объектов сплошной среды. Среди них - гибкий стержень, нить, контактный микропереключатель, пластина, мягкая оболочка, трехмерный упругий массив, явление цунами, вибрационная сепарация компонентов трехфазной смеси, гравитирующие жидкости в невесомости, сжатый газ под тяжелым поршнем и др.. [1-20]
Заключение
Для всех типов сплошной среды на основе физической дискретизации объекта исследования предложена обобщенная концепция формирования агрегированной математической модели его механического состояния относительно 6*N медленных переменных - координат собственных реперов у N дискретных элементов.
Предложенная научная технология обеспечивает:
• априорное разделение обобщенных координат на быстрые и медленные переменные;
• предварительные аналитические оценки инерционных и физических сил, входящих в ММ, по величинам медленных переменных;
• асимптотическое повышение точности результатов с увеличением числа элементов и информационную близость расчетной картины к реальному поведению объекта;
• физическую наглядность и низкую трудоемкость реализации математической модели;
• сведение технологии компьютерного анализа задач механики к задачам динамики системы абсолютно твердых тел;
• применимость алгоритмов численного анализа систем абсолютно твердых тел к задачам механики для сплошных сред.
Литература
1. Русанов П.Г. Построение расчетных моделей динамики сплошной среды с помощью метода твердых тел. -М.: МАИ: сб. докл. Научн. симпоз. "Динамические и технологические проблемы механики конструкций и сплошных сред", 1995, с. 41-42.
2. Русанов П.Г. Анализ динамики контактного микропереключателя методом твердых тел. -М.: МАИ: сб. докл. Научн. симпоз. "Динамические и технологические проблемы механики конструкций и сплошных сред", 1995, с. 40-41.
3. Русанов П.Г. Применение метода твердых тел в динамике деформируемого тела. Сб. докл. 10-той Зимней школы Ур. Отд. РАН по механике сплошных сред, Пермь: 1995, с. 213-214
4. Русанов П.Г. Моделирование цунами методом твердых тел. -Пермь: сб. докл. 10-той Зимней школы Ур. Отд. РАН по механике сплошных сред, 1995, с. 213
5. Русанов П.Г., Русанова Е.М. Компьютерный анализ динамики хлопка полукольцевой арки симметричном нагружении силой в плоскости кольца. Ш-межд. семинар" Технологические проблемы прочности", Подольск, 1996, с. 146-148.
6. Русанов П.Г., Русанова Е.М., Тихонова О.Н. Продольные колебания однородного вращающегося упругого стержня. Расчет статических смещений сечений и собственных частот колебаний методом твердых тел. -Подольск.: Госкомобр. 4-тыйМежд. сем. Технологические проблемы прочности, 1997, с 231-235.
7. Русанов Г.П., Гаврюшин С.С. Русанов П.Г. Численное моделирование динамики микропереключателя с упругими элементами. У-межд. сем. "Технологические проблемы прочности" Подольск, 1998, с. 193-200.
8. Жуков И.В., Кетат В.В., Русанов П.Г. Анализ динамики упругой балки методом физической дискретизации. Изв. Вузов. "Машиностроение". 2000. №3, с. 3-9.
9. Русанов П.Г. Анализ статики пластины, нагруженной неконсервативными силами, методом твердых тел. - Подольск.: Госкомобр. 8-ой Межд. сем. "Технологические проблемы прочности", 2001. с. 188-197.
10. Русанов П.Г. Динамика нити в однородном потоке воздуха. Рук. Х-межд. сем. "Технологические проблемы прочности", Подольск, 2003, с. 158-162.
11. Русанов П.Г. Отличия расчетов формы упругой линии гибкого стержня методами дискретизации МКЭ и МТТ. Изв. Вузов. "Машиностроение". 2006. № 8, с. 3-9.
12. Русанов П.Г. Численное моделирование динамики симметричного раскрытия мягкой оболочки в поле силы тяжести. Изв. Вузов. Машиностроение, 2006, № 10. с. 3-10.
13. Русанов П.Г. Динамика слияния двух капель в невесомости. ХЫУ Всеросс. Конф. по проблемам математики, информатики, физики и химии, РУДН, 2008. с. 37-39.
14. Русанов П.Г. Формирование, отрыв, полет и удар капли. ХЫУ Всеросс. Конф. по проблемам математики, информатики, физики и химии, РУДН, 2008. с. 39-40.
15. Русанов П.Г. Расчет собственных колебаний в вертикальной плоскости отрезка тяжелой нити. Изв. Вузов. Машиностроение, 2008, № 1, с. 3-10.
16. Русанов П.Г. Расчет собственных колебаний отрезка тяжелой нити с закрепленными концами «из вертикальной плоскости». Изв. Вузов. Машиностроение, 2008, № 2, с. 1-11.
17. Русанов П.Г. Влияние упругого звена на частоты собственных колебаний вертикально расположенного отрезка тяжелой нити. Изв. Вузов. Машиностроение, 2008, № 6, с. 3-7.
18. Русанов П.Г. Частоты собственных колебаний отрезка эластичной нити с неподвижным
концом. Изв. Вузов. Машиностроение, 2008, № 8, с. 3-9.
19. Русанов П.Г. Частоты собственных колебаний в вертикальной плоскости консольно закрепленного гибкого стержня. Изв. Вузов. Машиностроение, 2010, № 11, с. 3-7.
20. Русанов П.Г. Спектр частот собственных колебаний консольного гибкого стержня с вертикальной осью закрепленного участка. -М.: Вестник МГОУ, 2011, № 1, с. 17-24.
