Содержание к диссертации
Введение
ГЛАВА 1. Численный анализ контактных задач механики деформированного твердого тела 13
1.1. Математическая постановка контактной задачи теории упругости 13
1.2. Альтернирующий метод Шварца 15
1.3. Вариационная постановка задачи теории упругости 17
1.4. Вывод соотношений метода конечных элементов 19
1.5. Алгоритм численного решения контактной задачи теории упругости .39
1.6. Выбор итерационных параметров 48
1.7. Решение контактной задачи с учетом трения 54
1.8. Учет вязкоупругопластического деформирования 57
1.9. Решение уравнений МКЭ 63
1.10. Учет кинематических граничных условий при решении уравнений .66
1.11. Выводы к первой главе 70
ГЛАВА 2. Программная реализация алгоритма численного анализа контактных задач 71
2.1. Общая структура комплекса прикладных программ 71
2.2. Программы подготовки данных 73
2.2.1. Программы построения конечно-элементной модели 74
2.2.2. Программы аппроксимации физико-механических свойств материалов расчетной схемы 79
2.2.3. Программы формирования граничных условий 82
2.3. Программы центрального вычислительного блока 84
2.3.1. Основные процедуры при решении контактных задач теории упругости 84
2.3.2. Основные процедуры при решении вязкоупругопластических контактных задач 89
2.4. Программы представления данных 97
2.5. Выводы ко второй главе 102
ГЛАВА 3. Численные исследования контактного взаимодействия 103
3.1. Контактное взаимодействие полупространства и цилиндра 103
3.2. Контактное взаимодействие полупространства и полуцилиндра 106
3.3. Контактное взаимодействие неравномерно нагретых пластин 110
3.4. Поликонтактное взаимодействие 113
3.5. Напряженно-деформированное состояние резьбового соединения 117
3.6 Выводы к третьей главе 120
Основные результаты работы 121
Литература
- Вариационная постановка задачи теории упругости
- Учет кинематических граничных условий при решении уравнений
- Программы аппроксимации физико-механических свойств материалов расчетной схемы
- Контактное взаимодействие неравномерно нагретых пластин
Вариационная постановка задачи теории упругости
Многие ответственные узлы и элементы конструкций объектов энергетического оборудования, авиационной, аэрокосмической, наземной и морской техники работают в условиях контактного взаимодействия. Для оценки их ресурса и надежности требуется определять напряженно-деформированное состояние, для чего необходимо решать соответствующую контактную задачу. Вообще, контакт - это наиболее распространенный способ приложения нагрузок к деформируемому твердому телу, и концентрация напряжений в контактной зоне во многих случаях становится причиной разрушения материала. Как следствие, контактные задачи входят в число важнейших разделов в механике деформируемого твердого тела.
Основополагающие работы в теории контактного взаимодействия принадлежат Г. Р. Герцу, который получил распределение напряжений в зоне контакта упругих тел. Заметный вклад в развитие методов и алгоритмов аналитического решения контактных задач содержится в трудах советских ученых - И. Н. Векуа, Н. П. Векуа, Н. И. Мусхелишвили, Л. А. Галина, С. Г. Михлина, В. Л. Рвачева, Д. И. Шермана, И. Я. Штаермана, и многих других, а также работах зарубежных механиков и математиков К. Каттанео, Д. Синьорини, Р. Д. Миндлина, Н. Губера и других.
В связи с важностью и сложностью контактных задач, они до настоящего времени продолжают привлекать большое число исследователей в России (И.И. Ворович, В.М. Александров, В.И. Моссаковский, А.В. Манжиров, СМ. Айзикович, B.C. Давыдов, А.С. Кравчук, М.И. Чебаков, И.И. Аргатов, А.Г. Горшков, Н.Н. Дмитриев, Е.Н. Чумаченко, А.А. Успехов, Э.Р. Гольник и др.), и за рубежом (A. Curnier, G. Pietrzak, М. Barboteu, Е. Petoch, F. Armero, F. Lebon, E. Zahavi, D. Barlam, P. Wriggers и др.).
Для большого числа видов контактного взаимодействия и форм контактирующих поверхностей, в частности для большинства практически важных задач, не получено соответствующих аналитических решений. В этих случаях наиболее перспективными представляются численные методы [104, 105, 108, ПО, 111]. В задачах, характеризующихся сравнительно невысокими требованиями к гладкости входящих в их формулировку функций, основным методом расчетов является метод конечных элементов (МКЭ) [1 - 12].
Основными достоинствами МКЭ, способствующими его повсеместному применению, являются, например, высокая технологичность доступность и универсальность. МКЭ позволяет проводить численный анализ в областях сложной геометрической формы, учитывать граничные условия различных типов, а также физико-механические свойства материалов расчетных схем. Основные вычислительные процедуры МКЭ просты и прозрачны, что позволяет эффективно отслеживать обработку данных. Наконец, метод конечных элементов как алгоритмически, так и программно удобен для объединения с современными и средствами визуализации.
Повсеместное использование метода конечных элементов в процессе комплексной автоматизации сквозного цикла: проектирование -конструирование - изготовление, привело к возникновению многочисленных комплексов и пакетов прикладных программ, которые можно условно разделить на исследовательские и профессиональные. Для исследовательских КПП наибольшее значение имеет высокая скорость разработки, отладки и проведения численных исследований. В большинстве случаев, они сопровождаются непосредственно авторским коллективом и обладают узкоспециализированной ориентацией. После доработки пользовательского интерфейса, а иногда, и функционального ядра с учетом опыта эксплуатации, возможно доведение исследовательских КПП до профессионального уровня. Примерами известных профессиональных программных пакетов являются такие «средние» и «тяжелые» CAD/CAE/CAM системы, как отечественные МАРС, АСТРА, ЛИРА, FEMHCA, МЕГРЭ-ЗД, КАСКАД-2, МАК, ASTA, и зарубежные ANSYS, CATIA, MSC/NASTRAN, I-DEAS, MATRA, Pro/ENGINEER, Unigraphics, MARC+Mentat II, STAADIII, COSMOS/M, BEASY и др. Такие ППП благодаря своей высокой универсальности позволяют решать абсолютное большинство задач, не содержащих специфических сложностей. Основной трудностью в процессе их эксплуатации является достаточно высокая сложность овладения навыками сопровождения, причем процесс полного освоения профессионального пакета, как правило, оказывается длительным и трудоемким. Сравнительный анализ CAD/CAE/CAM систем приводится в работах [13-16].
В современных САЕ-системах, таких как ANSYS, LS-DYNA, MSC/NASTRAN, численное решение контактных задач производится, как правило, в рамках конечно-элементной технологии, при этом применяются следующие алгоритмы: метод внутренних многоточечных связей (МРС Algorithm); метод множителей Лагранжа (Pure Lagrange multiplier method); расширенный метод Лагранжа (Augmented Lagrange Method), метод штрафных функций (Penalty Method); комбинированный метод штрафов и Лагранжа (Lagrange&Penalty Method).
В число основных алгоритмов решения контактных задач входят метод штрафов, метод множителей Лагранжа и их комбинации [87 - 102], релаксационные схемы [17 - 20, 22], а также методы, основанные на введении контактных конечных элементов [98] или «псевдосреды» [21].
Из числа перспективных, но недостаточно разработанных методов решения контактных задач МДТТ выделяется альтернирующий метод Шварца, основанный на принципе поочередности [23 - 34]. Этот метод требует итерационного уточнения границ зон контакта, однако при его применении не требуется переформирование матриц систем линейных алгебраических уравнений. В отличие от большинства других методов, он также, как правило, не требует при построении конечно-элементных моделей согласовывать расположение узлов на контактных поверхностях.
В настоящее время наблюдается рост требований к точности и эффективности производимых численных исследований, а также продолжается расширение возможностей вычислительных средств. Поэтому сохраняется необходимость создания и дальнейшего развития методов решения контактных задач МДТТ, реализующих их алгоритмов и исследовательских КПП. Практика численных исследований убедительно показывает, что, одновременно с созданием новых и развитием существующих КПП общего назначения, требуется вести разработку целевых программ, направленных на решение задач в рамках одной или нескольких близких математических моделей, поскольку такие программы в своей предметной области позволяют существенно увеличить эффективность вычислительного эксперимента [36].
Применение метода конечных элементов к новым задачам требует существенного изменения отдельных этапов его реализации, в том числе перестроения собственно математических моделей, описывающих сложные физико-механические процессы, и создание новых численных алгоритмов их реализации для получения достоверных результатов. При этом организация вычислительных процедур должна проводиться с достаточно высокой точностью и максимальной экономичностью.
Актуальность проблемы. Интенсивное развитие методов математического моделирования как эффективного средства исследования сложных процессов деформирования с учетом контактного взаимодействия является одной из актуальных проблем прикладной математики, так как открывает новые возможности в развитии таких предметных областей как механика деформируемого твердого тела и прикладные методы численного анализа, значительно расширяет перспективы создания и практического использования систем автоматизированного проектирования.
С точки зрения прикладной математики особенно важным является дальнейшее развитие перспективных методов математического моделирования, применяемых для решения новых классов задач вычислительной термомеханики, математические постановки которых в общем виде учитывают сложные физико-механические эффекты, возникающие при вязкоупругопластическом деформировании с учетом контактного взаимодействия и сложного термосилового нагружения. Это дает возможность осуществления более полного анализа напряженно-деформированного состояния ответственных элементов конструкций, подверженных неоднородному термосиловому нагружению, и увеличения точности оценки их ресурса.
Цель работы. В соответствии с изложенным выше целью настоящей диссертационной работы является развитие перспективных численных методов решения квазистатических нелинейных краевых задач вычислительной термомеханики, учитывающих особенности контактного взаимодействия ответственных элементов конструкции в условиях сложного термосилового нагружения.
Учет кинематических граничных условий при решении уравнений
Последнее слагаемое в сумме (1.33), может быть учтено следующим образом. Аппроксимируем поверхность Sk, непосредственно участвующую в контакте, «поверхностными» одномерными конечными элементами аналогично поверхности S2. Тогда, с учетом соотношений (1.44) и (1.45), получим где Kk - число «поверхностных» конечных элементов, аппроксимирующих поверхность Sk. Вектор контактной нагрузки IRJf \ имеет вид компоненты контактной нагрузки, отнесенные соответственно к узлам грани / и j. Эти компоненты корректируются в процессе решения контактной задачи. Более предпочтительный подход состоит в следующем. Рассмотрим по отдельности каждый узел поверхности контакта Sk, непосредственно участвующих в контакте. Общее число таких узлов равно Кк. Член Ц., учитывающий контактное взаимодействие, в выражении (1.15) можно записать в виде где [i ] - глобальная матрица жесткости; {і?} - глобальный вектор нагрузки. Из этих соотношений видно, что глобальный вектор нагрузки {і?} и глобальная матрица жесткости [к] представляют собой суммы соответствующих локальных векторов и матриц с учетом геометрических связей, определяемых конечно-элементной моделью.
Таким образом, глобальная СЛАУ в матричной форме записывается в виде [K]{U} = {R}. (1.59) При решении системы (1.59) требуется учитывать кинематические граничные условия (1.4) и кинематические условия сопряжения на поверхности контакта Sk (1.6) или (1.9), если параметр /3 = 0. В случае решения системы прямым методом, для учета этих условий необходимо выполнить некоторые преобразования. В данной работе система уравнений (1.59) была решена при помощи итерационного (градиентного) метода и указанные условия учитывались в процессе решения.
Нужно отметить, что матрицы геометрических связей _ 2(е)_, _&(g)J и не используются в реальных вычислениях непосредственно. Глобальные матрицы, присутствующие в выражениях (1.52) - (1.57), формируются с помощью аппарата списков связности (матриц-указателей), построение и применение которых описано в работе [39]. 1.5. Алгоритм численного решения контактной задачи где р и q - компоненты вектора контактных сил в узле т (\ т Ма).
В соответствии с алгоритмом компоненты векторов \UЛ и \Rk}( ,, а є {А, В} поочередно и последовательно корректируются. Коррекция компонент вектора \Uk\ , ає{А,В} происходит на каждой четной итерации; коррекция компонент вектора \Rk\( ,, ссє{А,В} - на каждой нечетной. Рассмотрим произвольный контактный узел т (\ т МА , для определенности, лежащий на поверхности Sf тела А При отсутствии начального зазора коррекция компонент вектора \UЛ производится по формуле контактной поверхности Sf тела 5; ос 1т - итерационный параметр.
Выбор значений компонент и0 и v начального приближения вектора перемещений (и = 0) в случае контакта двух тел оказывает незначительное влияние на количество итераций. Значения выбираются на основании априорной информации о характере контактного взаимодействия.
Компоненты вектора контактных узловых сил \Rk\, (здесь и далее под узловыми силами подразумеваются поверхностные силы, приведенные к узлам при помощи функций формы конечного элемента) корректируются по формуле вектор узловых сил, приложенных в сходственной точке S контактной поверхности Sk тела В; сс(1т итерационный параметр. и вектора перемещении Uf = {U}(B) s = \ \ в сходственной точке s, лежащей на поверхности контакта Sk тела В. Построение сходственных точек может осуществляться разными способами [29 - 30, 34]. Например, для построения сходственной точки s можно опустить перпендикуляр из узла т на поверхность контакта «Sf тела В. Однако это, в общем случае, требует построения нескольких перпендикуляров к «поверхностным» конечным элементам, аппроксимирующим поверхность Sk и выбор основания, лежащего в пределах соответствующей грани. Вместо этого, строится отрезок т0т, соединяющий начальное т0 и конечное т положения узла, принадлежащие начальной Sk0 и конечной Sf контактным поверхностям тела А соответственно. Затем находится точка его пересечения с поверхностью тела В, лежащая на отрезке ij , где і и j - контактные узлы, принадлежащие поверхности тела В (Рис. 1.8).
Программы аппроксимации физико-механических свойств материалов расчетной схемы
Структурная схема программного комплекса представлена на Рис. 2.1. Программы можно разделить на три большие группы (блока): препроцессор, процессор и постпроцессор. Первый предназначен для подготовки данных, используемых при решении контактных задач МДТТ. Блок процессора составляют программы формирования глобальной СЛАУ, а также программы, предназначенные для решения СЛАУ с учетом кинематических граничных условий (в том числе контактных). Постпроцессор состоит из программ, используемых для представления результатов обработки данных при подготовке и решении контактных задач. Все программы комплекса, кроме специально оговоренных случаев, написаны на алгоритмическом языке ФОРТРАН 90. Программный код имеет объем около 25000 строк.
Обмен информацией между группами программ и программами внутри групп производится в виде записанных на винчестер рабочих файлов, образующих базу данных. Используемые объемы внешней и оперативной памяти согласовываются в соответствии с имеющимися вычислительными средствами.
Работа блоков КПП производится последовательно в соответствии алгоритмом решения задачи. Каждый блок вызывается в оперативную память процессора с внешнего носителя, после чего производится обработка данных и результаты заносятся в зону рабочих файлов в качестве исходной информации для следующих блоков. Построение комплекса таким образом позволяет, в случае необходимости, корректировать исходные данные и повторять некоторый шаг задания, полностью сохраняя все результаты работы предыдущего блока программ. В число приоритетов в процессе разработки алгоритмов и программ входило сокращение объемов исходных данных, так как при подготовке численных исследований это существенно уменьшает число ошибок, и увеличению универсальности вспомогательных алгоритмов, предназначенных для наиболее полного учета различных сложных кинематических и силовых граничных условий.
Разработанный программный комплекс является исследовательским, однако он обладает достаточными возможностями в области 2D-моделирования, чтобы послужить основой для создания «легкой» САЕ-системы. Комплекс применялся для проведения численных исследований ответственных элементов конструкций, находящихся в состоянии контактного взаимодействия [69, 72]. Результаты некоторых из них отражены в третьей главе данной работы.
В группу препроцессора входят программы, предназначенные для построения сетки конечно-элементной модели, массовых сил, кинематических и силовых граничных условий (поверхностной и узловой механической нагрузки), физико-механических свойств расчетной схемы. К физико-механическим свойствам относятся компоненты тензора линейного расширения, модули упругости, коэффициенты Пуассона, данные для аппроксимации диаграммы деформирования и т. п.
Поверхностная нагрузка, в общем случае неравномерная, задается на сторонах конечных элементов («поверхностных» элементах), аппроксимирующих внешнюю поверхность исследуемых тел (Таблица 3). Узловая нагрузка рассматривается в виде дискретных сил, которые прикладываются строго в узлах «поверхностных» конечных элементов (Таблица 4). Силы тяжести задаются в качестве массовых сил. Также в некоторых расчетных схемах присутствуют неравномерные тепловые поля, значения которых задаются в узлах или на элементах.
Данные группируются, аппроксимируются и записываются в зависимости от характера и объема вводимой информации в разные области зоны рабочих файлов, расположенной на внешнем носителе. Предусмотрен табличный вывод для контроля вводимых данных и графическая визуализация некоторых параметров.
В качестве математического обеспечения при подготовке данных используются сплайн-аппроксимация или полиномиальная аппроксимация, основанная на методе наименьших квадратов.
Программы построения конечно-элементной модели Автоматической генерации конечно-элементных сеток в двухмерных областях при помощи трехузловых симплекс-элементов посвящено большое количество работ, в том числе [73, 74, 76]. В данной работе для численного решения контактных задач использовался алгоритм построения двухмерных конечно-элементных моделей, описанный в данном разделе [75]. Он реализуется следующим образом. Первым шагом является создание геометрической двухмерной модели исследуемой конструкции с помощью «легкой» CAD-системы типа AutoCAD или CADKEY. Использование более сложных систем, например SolidWorks, КОМПАС-ЗБ, T-FLEX CAD, SolidEdge, также позволяет получать требуемые изображения, однако представляется излишним, так как полученную геометрическую модель для триангуляции требуется перестроить более удобным образом с учетом нагружения. В данной работе был использован пакет AutoCAD 2010.
Следующим шагом после адаптации является разбиение геометрической модели на конечные элементы. В «средних» и «тяжелых» САЕ-системах для этого предусмотрены специальные ППП, например, PATRAN и FEMAP, имеющие генераторы двухмерных и трехмерных конечно-элементных сеток, в том числе в областях достаточно сложной геометрической формы.
В данной работе реализованы алгоритмы и программы, предназначенные для автоматического построения сеток из трехузловых симплекс-элементов в двухмерных областях произвольной геометрической формы, в том числе многосвязных [75], на основе данных о внутренних и внешних контурах разбиваемых областей. Исходная область разбивается на подобласти, в нескольких точках контуров которых задается допустимый размер шага по пространству. Для каждой области определяются длины контуров и между заданными точками с учетом требований к размеру шагов размещаются контурные узлы. После разбиения всех контуров подобласти последовательно триангулируются в соответствии со следующим алгоритмом. Рассмотрим произвольный контур, содержащий некоторое число контурных узлов (Рис. 2.2).
Контактное взаимодействие неравномерно нагретых пластин
Для иллюстрации работы алгоритма расчета поликонтактного взаимодействия в условиях вязкоупругопластического деформирования было исследовано контактное взаимодействие четырех пластин, подвергнутых неравномерному температурному нагружению (Рис. 3.17).
Расчетная схема поликонтактного взаимодействия Материал пластин - сталь 40Х. Размер каждой из пластин составляет 0.03x0.025 м. На верхнюю пластину действует распределенная вертикальная нагрузка р2 = -200 МПа. Температура в пределах всех пластин меняется линейно вдоль оси Охх от оси закрепления Ох2 до правой свободной поверхности (0 хх 0.03). В точках обоих тел, лежащих на оси Ох2, задана температура 600К, а в точках свободной поверхности ( = 0.03) -температура 200К.
Полученная геометрия конечно-элементных моделей пластин после деформирования приведена на Рис. 3.18 (перемещения увеличены). Видно, что характер деформирования определяется в первую очередь температурным воздействием. Распределение компоненты т22 тензора напряжений показано на Рис. 3.19. Найденное поле напряжений имеет неоднородный характер и в зонах непосредственного контакта аналогично полученному в задаче Герца. Максимальные сжимающие напряжения, порядка 900 МПа, наблюдаются на контактных поверхностях, что приводит к возникновению незначительной пластической деформации. Напряжения в зонах растяжения составляют не более 90МПа. На Рис. 3.20 представлено поле компоненты є22 тензора полной деформации. Поле компоненты т99 Рис. 3.20. Поле компоненты є тензора напряжений, МПа тензора деформации, 10
Также исследовано изменение напряженно-деформированного состояния пластин под действием ползучести. На Рис. 3.21, 3.23, 3.25 показано распределение компоненты т22 тензора напряжений в разные моменты времени. На Рис 3.22, 3.24, 3.26 показано распределение компоненты є22 тензора полной деформации. Время, прошедшее с начала процесса равно 100 ч, 200 ч и ґ = 400ч соответственно. Видно, что пластины подвергаются одновременно эффектам ползучести и релаксации.
Напряженно-деформированное состояние резьбового соединения Численное решение осесимметричной контактной задачи МДТТ проиллюстрировано на примере резьбового соединения (Рис. 3.27). Материал конструкции - сталь 40Х. Геометрические размеры соединения взяты из работы [86]. Предполагалось, что растягивающая нагрузка, приложенная к конструкции, вызывает известное фиксированное перемещение и2 = 1 10 6 м.
Разработанный численный алгоритм и реализующий его КПП протестированы на задаче Герца. Полученные результаты хорошо согласуются с известным аналитическим решением.
Рассмотрены плоские и осесимметричные краевые контактные задачи МДТТ в областях сложной геометрической формы, демонстрирующие возможности разработанного алгоритма. Получены результаты расчетов с учетом поликонтактного взаимодействия и термовязкоупругопластического деформирования.
1. На основе альтернирующего метода Шварца разработаны математические модели и итерационные алгоритмы численного исследования квазистатических краевых контактных задач механики деформируемого твердого тела в двухмерных областях, имеющих сложное геометрическое оформление, с учетом вязкоупругопластического деформирования и неоднородного термосилового нагружения.
2. С применением разработанных алгоритмов создан комплекс прикладных программ, предназначенный для решения контактных задач МДТТ в двухмерных областях.
3. Получены решения некоторых двухмерных контактных задач механики деформируемого твердого тела, демонстрирующие возможности разработанных математических моделей, алгоритмов и программ.
