Numerical Simulation of Large Plastic Deformations and Free Flight of a Beam After Impact with a Rigid Stop

Full Text

1. Введение Современные численные методы позволяют решать все более сложные статические и динамические задачи [1-5], при этом учитываются многие факторы, пусть и усложняющие расчет, но делающие расчетные модели более общими, приближающими результаты анализа к поведению реальных объектов. Данная работа использует метод конечных элементов и явную вычислительную схему интегрирования уравнений движения, что в настоящий момент не является безусловным признаком новизны и ранее исследовалась в [6-10]. Вместе с тем в настоящей работе авторами применяются только собственные вычислительные алгоритмы [11-15] и собственная компьютерная программа, что дает определенные преимущества среди которых полное знание деталей алгоритма, возможность легкого внесения необходимого изменения в программные модули, иногда ориентированные даже на конкретные задачи. Использование собственной компьютерной программы, конечно, не гарантирует научную новизну, многие элементы алгоритмов повторяют общепринятые подходы. Наибольшее внимание уделено тем элементам компьютерной программы и алгоритма, которые авторы считают новыми: работа с напряжениями в глобальной системе координат, несмотря на рассмотрение больших перемещений, и углов поворота; контроль результатов решения с использованием понятия «центра масс»; использование равномерной конечно-элементной сетки в начале решения, возможность рассмотрения пластических сдвигов лишь на тех площадках, где превышен предел упругости. В исследовании авторы ограничились рассмотрением плоской деформации, условием наступления текучести по максимальным касательным напряжениям. Предлагаемый алгоритм ориентирован на адекватное отображение следующих особенностей поведения расчетной модели: большие перемещения, углы поворота; большие деформации; отображение волновых процессов; отображение процессов разгрузки, если таковая имеет место в какой-либо области расчетной схемы; моделирование свободного полета; формирование контактных усилий с использованием закона трения по Кулону. Дальнейшее изложение будем вести на примере достаточно сложной модельной задачи, делая при необходимости теоретические отступления. Авторами используется простейший конечный элемент с линейным законом для перемещений и постоянными напряжениями внутри элемента. Недостатки этой модели хорошо известны и преодолеваются прежде всего в сгущении сетки. 2. Метод решения 2.1. Интегрирование уравнений движения Используется явная вычислительная схема интегрирования уравнений движения с экстраполяцией по Адамсу. Распределенная масса заменяется точечными массами в узлах конечно-элементной сетки. На каждом временном шаге необходимо лишь определить сосредоточенные силы, действующие на точечные массы. Применение 2-го закона Ньютона дает соответствующие ускорения. Дальнейшая процедура вычислений описана во многих наших публикациях [11-15], поэтому подробнее на ней останавливаться не будем. Явная схема не требует решения системы линейных уравнений, но требует использования крайне малых шагов по времени, что, как ни странно, можно отнести как к недостаткам, так и к преимуществам метода. Преимущества - гарантированное отражение волновых процессов и точная фиксация по времени перехода в пластическое состояние или выхода из него. 2.2. Описание расчетной схемы численного примера На рис. 1 представлена плоская расчетная схема. Приведем полные характеристики расчетной схемы. Материал балки (точнее балки стенки) максимально приближен к обычной стали 3, модуль упругости Е = 2 000 000 кг/см2, коэффициент Пуассона , удельный вес 0,00785 кг/см3, предел текучести при растяжении 2400 кг/см2. Коэффициент трения Кулона по контактной поверхности 0,35. Упор абсолютно жесткий, неподвижный. Скорость падения балки принята 80 м/с, что примерно соответствует установившейся скорости падения небольшого предмета в воздушной среде. До начала касания напряжения в пластине практически отсутствуют, но действует собственный вес. Внутреннее трение в материале не учитывается, если не считать, конечно, колоссальной потери энергии при пластической деформации, которая и гасит большую часть начальной кинетической энергии пластины. Рис. 1. Расчетная схема с самой грубой сеткой (800 конечных элементов) И с то ч н и к: выполнено А.В. Штейном Figure 1. Calculation model with the coarsest grid (800 finite elements) S o u r c e: compiled by A.V. Shtein 2.3. Моделирование контактных усилий Контактное взаимодействие может быть рассмотрено во многих, в том числе широко употребляемых в настоящее время программных конечно-элементных комплексах, например LS-DYNA, NASTRAN и др. В наших ранних публикациях [12; 15] приводится вариант вычисления контактных усилий, который используется в данной работе. В последней из упомянутых публикаций контактные силы вычисляются для трехмерного случая. Чтобы не повторять содержание ранее опубликованного материала, ограничимся лишь кратким изложением учета контактного взаимодействия. Касание внешних узлов одного тела (в нашем случае это будет пластина) с участками границы другого тела (в нашем случае это будет упор) осуществляется через некоторые пружины (контактные жесткости). Нормальная составляющая контактной силы считается пропорциональной глубине проникновения внешнего узла вглубь границы тела. Касательная составляющая контактной силы накапливается при смещении узла вдоль линии контакта. Когда она достигает предельного значения по Кулону, начинается скольжение по линии соприкосновения. Контактное взаимодействие может наступать, прекращаться и снова возникать в процессе счета. Все эти операции естественно выполняются автоматически. На вычисления контактных сил тратится определенное процессорное время. С увеличением жесткости контактных пружин точность моделирования контактного взаимодействия повышается, так как уменьшается зазор между соприкасающимися телами, который в идеале должен вообще отсутствовать. Однако увеличение контактных жесткостей приводит к необходимости уменьшения шага интегрирования, поэтому здесь приходится идти на некоторый компромисс. Касательная и нормальная жесткости для рассматриваемых примеров были приняты равными 2 × 106 кг/см и для всех типов сетки. Поскольку с уменьшением размера конечного элемента количество пружин на единицу длины контактного взаимодействия будет возрастать, можно говорить о том, что со сгущением сетки конечных элементов точность контактного взаимодействия также возрастает. 2.4. Порядок вычисления узловых реакций в конечном элементе Использован простейший треугольный конечный элемент с линейным полем перемещений и соответственно постоянными по объему напряжениями. Когда напряжения в элементе известны, узловые реакции легко определяются. На этой стандартной для метода конечных элементов процедуре не будем останавливаться, сосредоточившись на алгоритме вычисления напряжений в элементе. В нашей модели напряжения не могут быть получены непосредственно по значениям деформаций, но зависят от истории их накопления. Будем далее исходить из того, что шаги по времени весьма малы, соответственно, малыми являются и приращения деформаций на каждом шаге. Сначала по приращениям смещений узлов на данном шаге определяются соответствующие приращения напряжений. Это делается точно так же, как в обычной линейной задаче. Правда, здесь следует использовать координаты узлов, которые имеет система в рассматриваемый момент времени. Полученные приращения напряжений добавляются к напряжениям текущего шага и получаются исправленные напряжения. Здесь удобно, чтобы напряжения и их добавки были записаны для одной системы координат, что и реализовано в нашем алгоритме. Далее вычисляются главные напряжения и выполняется проверка на превышение максимальных касательных напряжений: ; (1) ; (2) , (3) где , , - значения главных напряжений текущего шага, - предел текучести материала. Как было сказано выше, в данной работе авторы ограничились критерием пластичности наибольшего касательного напряжения. Если все три неравенства выполняются, это будет означать, что на данном шаге конечный элемент работает как упругий и не нужно вносить корректировку в значения напряжений. Если неравенство (1) не выполняется (практически наиболее частый случай), следует уменьшить и на столько же увеличить , чтобы неравенство (1) превратилось в равенство. Далее необходимо проверить остальные неравенства и в случае их нарушения (весьма редкий случай) соответствующим образом откорректировать значения главных напряжений. Если нарушение неравенств имело место, следует зафиксировать, что на данном шаге в данном элементе происходят пластические деформации. После корректировки главных напряжений необходимо сделать очевидную корректировку напряжений на исходных площадках. На каждом шаге необходимо еще определить поворот элемента как жесткого целого по формуле , где u и v - приращения перемещений точек элемента вдоль горизонтальной и вертикальной оси на текущем шаге интегрирования. Над только что полученными значениями напряжений при помощи угла следует выполнить преобразование поворота. Таким образом, наш алгоритм настроен на получение напряжений , , в неподвижных осях координат x и y. После того как напряжения в элементе определились, можно найти узловые реакции в конечном элементе по известной методике, которую не будем повторять. 3. Результаты В разделе 3 подробно описан пример для численного анализа. Моделирование выполнялось для трех конечно-элементных сеток. Грубая сетка, размер катета треугольного элемента - 1 см, количество элементов - 800. Средняя сетка, размер катета треугольного элемента - 0,25 см, количество элементов - 12 800. Мелкая сетка, размер катета треугольного элемента - 0,0625 см, количество элементов - 248 000. Обозначение сеток грубая, средняя и мелкая будет далее использоваться в тексте. Сгущая сетку, авторы каждый раз уменьшали линейный размер элемента в 4 раза. Если говорить о времени, которое затрачивает компьютер на решение примера, то переход от грубой сетки к мелкой требует увеличения этого времени примерно в 4 096 раз. Здесь надо иметь в виду, что при переходе к более мелкой сетке не только возрастает количество элементов, но и уменьшается шаг интегрирования в 4 раза. Если при использовании грубой сетки измеряется секундами, то при использовании мелкой сетки требуются уже сутки работы компьютера. Эти данные соответствуют использованию обычного компьютера с тактовой частотой 3,6 GHz. Большое время тратится также на создание компьютерной анимации и графиков, если эти действия выполняются. 3.1. Сходимость при сгущении сетки конечных элементов Сначала рассмотрим сходимость по пространственной координате. В табл. 1 показаны состояния системы в момент времени t = 0,0005 c момента касания упора. Как можно видеть, три сетки дают близкий результат с точки зрения общего состояния балки. Удовлетворительно совпадают области касания. Светлым цветом на рисунках отмечены элементы, которые в текущий момент времени находятся в состоянии пластического течения. Области, захваченные пластическими деформациями для средней и мелкой сетки, удовлетворительно совпадают. Таким образом, данные табл. 1 свидетельствуют о наличии сходимости по пространственной координате. В таблице даны значения шагов по времени, которые использованы при численном интегрировании. К сожалению, в статье нет возможности воспроизвести компьютерную анимацию, которая реализована в программе и дает лучшее представление о характере исследуемого движения, чем отдельные рисунки. В табл. 2 представлена балка, находящаяся в состоянии свободного полета после отскока от упора. Как видно, очертания графиков (исключая ускорения, их обсудим позже) визуально совпадают. На каждом графике имеются три цифры - максимальное, текущее и минимальное значение, что позволяет установить масштаб графика. Контроль точности удобно производить по максимальным и минимальным значениям. Максимальные перемещения для двух сеток отличаются на 8 %, значения скорости - на 2 %. Сходимость для ускорений существенно хуже, максимальные значения отличаются на 34 %, очертания графиков далеко не полностью совпадают. Вообще при сгущении сетки ускорения в конкретной точке при ударных воздействиях неограниченно возрастают. Метод численного интегрирования, который здесь используется, отражает и самые высокочастотные колебания, которые происходят в расчетной схеме. Чем меньше размеры конечных элементов, тем больше будут частоты соответствующих высших форм и больше ускорение конкретного узла расчетной схемы. Эта особенность рассмотрена в нашей публикации [14]. Таким образом, использование ускорения узла для анализа сходимости не имеет в данном случае смысла. Таблица 1 / Table 1 Область, захваченная пластическими деформациями в момент времени 0,0005 с с начала касания упора при различных размерах конечных элементов / The area captured by plastic deformations at a time of 0.0005 s from the beginning of the contact with the stop at various sizes of the finite elements Название сетки Кол-во конечных элементов / Name of the grid Number of finite elements Состояние пластины и область, где происходят пластические деформации в текущий момент времени / The state of the plate and the area where plastic deformations occur at the current time Грубая сетка 800 элементов / Coarse grid of 800 elements Шаг по времени Dt = 0,000000256 с / Time step Dt = 0.000000256 s Средняя сетка 12 800 элементов / The Medium grid of 12,800 elements Шаг по времени Dt = 0,000000064 с / Time step Dt = 0.000000064 s Мелкая сетка 204 800 элементов / Fine grid of 204800 elements Шаг по времени Dt = 0,000000016 с / Time step Dt = 0.000000016 s И с то ч н и к: выполнено А.В. Штейном S o u r c e: compiled by A.V. Shtein Таблица 2 / Table 2 Балка в состоянии свободного полета после отскока, момент времени 0,02 с c начала счета при различных размерах конечных элементов The beam in a state of free flight after rebound, at a time of 0.02 s since the beginning of counting for various sizes of finite elements Название сетки Кол-во конечных элементов / Name of the grid Number of finite elements Состояние пластины / The state of the plate Грубая сетка 800 элементов / Coarse grid of 800 elements Шаг по времени Dt = 0,000000256 с / Time step Dt = 0.000000256 s Средняя сетка 12 800 элементов / Medium grid of 12,800 items Шаг по времени Dt = 0,000000064 с / Time step Dt = 0.000000064 s Мелкая сетка 204 800 элементов / Fine grid of 204,800 elements Шаг по времени Dt = 0,000000016 с / Time step Dt = 0.000000016 s И с т о ч н и к: выполнено А.В. Штейном / S o u r c e: compiled by A.V. Shtein Для того чтобы все же дать оценку ускорений с точки зрения сходимости, можно предложить осуществлять слежение за ускорением центра масс балки. График ускорения центра масс балки приведен в нижней строке табл. 3. Как видно, здесь сходимость наблюдается. Вообще говоря, из расчетной схемы можно вырезать и некоторую область и рассматривать центр масс этой области. По существу, здесь получим некоторое усреднение ускорения по выбранной области. Далее рассмотрим сходимость для напряжений при сгущении сетки. Приведем графики первого главного растягивающего напряжения (табл. 4) для средней и мелкой сетки. Элементы, в которых определяются напряжения, показаны в нижней строке табл. 4. Таблица 3 / Table 3 Перемещение, скорость, ускорение левой верхней угловой точки пластины по вертикали, ускорение центра масс. Средняя и мелкая сетки / Vertical displacement, velocity, acceleration of the upper left corner point of the plate vertically, acceleration of the center of mass. Medium and fine grids Средняя сетка / Medium grid Мелкая сетка / Fine grid И с то ч н и к: выполнено А.В. Штейном / S o u r c e: compiled by A.V. Shtein Таблица 4 / Table 4 Первые главные напряжения в трех точках пластины / The first principal stresses at three points of the plate Средняя сетка / Medium grid Мелкая сетка / Fine grid И с т о ч н и к: выполнено А.В. Штейном / S o u r c e: compiled by A.V. Shtein Из анализа данных табл. 4 видно, что и для напряжений имеется сходимость при сгущении сетки. Существенно лучший результат получается для точки, лежащей на середине балки по высоте (расхождение менее 2 %), хуже результат для верхней и нижней точек. Объяснение здесь заключается в том, что верх и низ пластины - это области с резким изменением напряжений и здесь в полной мере проявляются дефекты треугольного элемента. Также надо иметь в виду, что фактически идет сравнение напряжений в разных точках. Различное положение точек слежения за напряжениями проиллюстрировано на рис. 2 для средней и мелкой сеток. Рис. 2. Положение граничной точки для средней и мелкой сетки И с т о ч н и к: выполнено А.В. Штейном Figure 2. The position of the boundary point for the medium and fine grid S o u r c e: compiled by A.V. Shtein Подводя итог, можно констатировать, что для напряжений сходимость, пусть не в лучшем виде, но достигнута. 3.2. Исследование сходимости при уменьшении шага интегрирования Далее перейдем к исследованию сходимости по временному шагу. Будем использовать тот же временной отрезок, что и при исследовании сходимости по напряжениям и среднюю сетку. В табл. 5 приведены численные результаты решения, полученные при трех значениях шага интегрирования , и . Приводятся максимальные по модулю значения, взятые с графиков. Кинематические параметры определялись для левой верхней точки пластины, точки для определения напряжений показаны в нижней строке табл. 4. Таблица 5 / Table 5 Исследование сходимости при уменьшении шага интегрирования по времени. Смещение, скорость и ускорение левой верхней точки пластины. Главные напряжения в трех точках / Investigation of convergence with decreasing integration step in time. Displacement, velocity, and acceleration of the upper left point of the plate. Principal stresses are at three points Шаг интегрирования, секунды / Integration step, seconds 6,410 ∙ 10 - 8 3,205 ∙ 10 - 8 1,300 ∙ 10 - 8 Смещение, см / Displacement, cm - 0,303 ∙ 10 - 0,303 ∙ 10 - 0,303 ∙ 10 Скорость, см/c / Velocity, cm/c - 0,906 ∙ 104 - 0,906 ∙ 104 - 0,906 ∙ 104 Ускорение, см/c2 / Acceleration, cm/c2 - 0,352 ∙ 109 - 0,348 ∙ 109 - 0,349 ∙ 109 верх, кг/см2 / upper, kg/сm2 -0,285 ∙ 104 -0,283 ∙ 104 -0,283 ∙ 104 середина, кг/см2 / middle, kg/сm2 - 0,106 ∙ 105 - 0,106 ∙ 105 - 0,106 ∙ 105 низ, кг/см2 / bottom, kg/сm2 - 0,360 ∙ 105 - 0,351 ∙ 105 - 0,354 ∙ 105 верх, кг/см2 / upper, kg/сm2 - 0,116 ∙ 104 - 0,119 ∙ 104 - 0,120 ∙ 104 середина, кг/см2 / middle, kg/сm2 - 0,130 ∙ 105 - 0,130 ∙ 105 - 0,130 ∙ 105 низ, кг/см2 / bottom, kg/сm2 - 0,384 ∙ 105 - 0,376 ∙ 105 - 0,378∙ 105 *Приведены максимальные по модулю значения за отрезок времени 0,00052 с / *The maximum modulo values for a time interval of 0.00052 s are given И с то ч н и к: выполнено А.В. Штейном / S o u r c e: compiled by A.V. Shtein Как видно с уменьшением шага по времени, расхождения даже для ускорения не превышают 0,3 %. Для напряжений расхождение составляет менее 0,8 %. Такую сходимость можно признать хорошей, значительно лучшей, чем по пространственной координате. 3.3. Анализ результатов с точки зрения моделирования волн деформаций В табл. 6 представлены деформированные состояния системы с эпюрой главного растягивающего напряжения . В статье приводятся рисунки только для нескольких моментов времени. Таблица 6 / Table 6 Демонстрация волнового характера решения, эпюры первого главного напряжения / Demonstration of the wave nature of the solution, a plot of the first principal stress Момент времени, с / Point in time, s Состояние системы / System state Шкала для напряжений / Scale for stresses t = 0,000010096 t = 0,000014135 t = 0,000018173 Продолжение табл. 6 / Continuation of the Table 6 Момент времени, с / Point in time, s Состояние системы / System state t = 0,000022212 t = 0,000036346 t = 0,000040365 t = 0,000044423 t = 0,000048462 Окончание табл. 6 / Ending of the Table 6 Момент времени, с / Point in time, s Состояние системы / System state t = 0,000405869 t = 0,000409908 И с то ч н и к: выполнено А.В. Штейном / S o u r c e: compiled by A.V. Shtein Первые четыре рисунка табл. 6 представляют начало процесса соударения. Используя их, можно получить скорость распространения фронта волны деформаций, которая оказывается близкой к аналитическому значению для неограниченного упругого тела: Необходимо подчеркнуть, что значение скорости распространения волны заметно не изменилось, несмотря на сильные отступления от закона Гука в модели. Также отметим, что в области касания возникает напряженное состояние, близкое к трехосному гидростатическому сжатию. Следующие четыре рисунка табл. 6 соответствуют известному правилу отражения волны сжатия от свободной поверхности волной растяжения. Последние два рисунка отражают стадию соприкосновения, за которой следует отскок. Сравнение с упругим решением показывает, что пластические деформации существенно удлиняют время касания при ударе, которое оказалось равно 0,00133 с. Если отказаться от учета пластических деформаций, время соприкосновения будет 0,00034 с. Соответствующий расчет авторами выполнялся, но в статье не приведен. 4. Заключение 1. Представлены основные этапы и последовательность решения динамической задачи о соударении упругопластических тел. Для решения использовалась компьютерная программа, полностью разработанная авторами статьи. Задача решена с учетом больших перемещений и деформаций, допускает рассмотрение свободного полета и контактного взаимодействия. Приведенное решение носит явно выраженный волновой характер. 2. Отличительные черты развиваемого подхода - получение и работа с напряжениями в неподвижной глобальной системе координат и возможность рассмотрения пластических сдвигов лишь по тем площадкам, где достигнут предел текучести, т.е. рассматривается сложное загружение. 3. Приведен пример численного моделирования для плоской задачи соударения пластины с упором, на основании которого удалось продемонстрировать удовлетворительную сходимость результатов при сгущении сетки конечных элементов и хорошую сходимость результатов при уменьшении шага интегрирования уравнений движения по времени. 4. Показано, что для используемого метода сходимость для ускорения точки при сгущении сетки не может быть достигнута и это не является признаком ущербности метода. 5. Предложено использовать понятие центра масс при анализе сходимости решения для конечно-элементных схем большой размерности.

About the authors

Alexey V. Shtein

Russian University of Transport

Email: avsh7714@yandex.ru
ORCID iD: 0009-0003-2232-5121
SPIN-code: 3150-4438

Associate Professor of the Department of Structural Mechanics

9/9, Obraztsova St, Moscow, 127994, Russian Federation

Vladimir B. Zylev

Russian University of Transport

Author for correspondence.
Email: zylevvb@ya.ru
ORCID iD: 0000-0001-5160-0389
SPIN-code: 5833-9120

Doctor of Technical Sciences, Professor, Department of Structural Mechanics

9/9, Obraztsova St, Moscow, 127994, Russian Federation

References

Supplementary files