Численное моделирование разделительных процессов обработки давлением Текст научной статьи по специальности «Математика»

УДК 621.96+539.3+519.6

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ РАЗДЕЛИТЕЛЬНЫХ ПРОЦЕССОВ

ОБРАБОТКИ ДАВЛЕНИЕМ

© 2012 К. М. Иванов1, П. М. Винник1, В. Н. Иванов2

балтийский государственный технический университет «Военмех»

имени Д. Ф. Устинова Санкт-Петербургский институт машиностроения (ЛМЗ-ВТУЗ)

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

Метод конечных элементов, разделительные операции, матрицы жёсткости, число обусловленности.

Введение

Решение современных задач технологической подготовки компьютерно-интегрированного производства требует наличия адекватных математических моделей. К числу наименее изученных с точки зрения численного моделирования процессов обработки давлением относятся разделительные операции.

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

Эти два обстоятельства затрудняют численное моделирование процесса даже с использованием самых современных программно-аппаратных средств.

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

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

Одним из основных методов численного моделирования процессов обра-

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

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

Причём по самой физической природе разделительных процессов предметом исследования являются расчётные поля тех или иных величин (деформаций, интенсивности деформации, напряжений и т.д.) именно в областях, в которых КЭ

претерпевают наиболее значительные изменения.

Возможны три подхода к этой проблеме.

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

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

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

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

В-третьих, таким образом можно, учитывая искажение формы некоторых КЭ, оценить влияние этого искажения на полученные расчётные поля. Для этого необходимо иметь возможность отслежи-

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

Нерегуляризованные матрицы жёсткости

Отметим, что искажение формы КЭ приводит к увеличению числа обусловленности локальной матрицы жёсткости элемента или, говоря другими словами, к тому, что некоторые собственные числа у этой матрицы становятся достаточно малыми по абсолютной величине. Но до регуляризации матрицы жёсткости часть собственных чисел локальной матрицы и так нулевые (например, в случае упруго-пластического анализа для случая плоской деформации локальная матрица жёсткости линейного треугольного КЭ имеет три нулевых и три положительных собственных числа). Таким образом, если для, например, задач без сильного искажения формы КЭ, в случае ещё не регуляризованной глобальной матрицы жёсткости [КО] нахождение перемещений узлов [и] приводит к системе линейных алгебраических уравнений (СЛАУ)

[КО] [и] = [ГС] (1)

(здесь [ТС] — вектор нагрузок), которая заданием некоторых перемещений приводится к СЛАУ

[К] [ио] = [Н (2)

с уже невырожденной матрицей жёсткости [К], то для рассматриваемых разделительных процессов в силу указанной возможной близости локальной матрицы того или иного элемента к вырожденной необходимо считаться с возможностью того, что даже после задания некоторых перемещений глобальная матрица жёсткости [К] может быть близка к вырожден-

ной, то есть может быть плохо обусловленной.

Таким образом, задача нахождения перемещений узлов метода конечных элементов при моделировании разделительных операций является, вообще говоря, некорректно поставленной задачей [1].

Заметим, что на простейших примерах легко показать, что выбор того или иного набора перемещений, задаваемых для регуляризации матрицы [КО], даже в случае получения невырожденной матрицы [К] изменяет её собственные числа (как следствие, и число обусловленности), то есть переход от матрицы [КО] к матрице [К] существенно неоднозначен. Так что с целью устранения этой неоднозначности целесообразно даже в этом случае иметь возможность оценить обусловленность по самой матрице [КО].

И уж для моделирования разделительных операций, когда есть все вышеуказанные основания ожидать получения вырожденной матрицы [К], тем более целесообразно не переходить к решению СЛАУ (2), а рассматривать непосредственно СЛАУ (1).

Как известно [1], решение СЛАУ с возможно вырожденной матрицей определено с точностью до вектора [¿] из множества Ы[Кс\ = {И: [К0]-[г]=[0]}, и поэтому полезно сначала оценить структуру множества Щко].

Так как глобальная матрица жёсткости [КО] является суммой (по всем КЭ) матриц, каждая из которых построена по локальной матрице жёсткости данного КЭ, то анализ структуры множества Щко] имеет смысл начинать с анализа структуры аналогичного множества для произвольного КЭ.

Можно показать, что для локальной матрицы жёсткости КЭ в рассматриваемом случае (упругопластический анализ для плоской деформации) трём нулевым собственным числам соответствуют три перпендикулярных друг другу собственных вектора (коэффициенты Ьи Ьр, Ьт, с, Ср ст понимаются традиционно для МКЭ),

Г1 > Г 0 ^ (Ъ. - ь Л } т

0 1 с. - с .т

1 0 Ъ - Ъ. т 1

0 1 с - с. т.

1 0 Ъ - Ъ.

V 0 0 V1 0 V с - с.,

имеющих следующий физический смысл (напомним, что какое-либо закрепление узлов КЭ отсутствует, так как не предполагается выполненной регуляризация матриц жёсткости): возможно произвольное смещение КЭ как твёрдого тела в направлении оси OX (первый вектор), оси OY (второй вектор) и произвольный поворот КЭ в плоскости XOY вокруг его центра тяжести (точки пересечения медиан) (третий вектор).

На этом основании можно показать, что для матрицы [КО] тремя базисными векторами множества NКО] будут векторы, имеющие тот же физический смысл, относящийся теперь уже ко всему множеству конечных элементов. Первые два вектора будут иметь аналогичную локальной матрице жёсткости структуру, а третий при более сложной структуре будет соответствовать повороту в плоскости XOY всего множества КЭ как одного твёрдого тела вокруг его центра тяжести.

Понятно, что напряжённо-деформированное состояние заготовки не зависит от таких перемещений заготовки как твёрдого тела.

Как известно, для положительно-определённой матрицы [А] спектральное число обусловленности с([А]) определяется как отношение максимального собственного числа к минимальному. Определим для неотрицательно-определённой матрицы [КО] спектральное число обусловленности с0([КО]) как отношение её максимального собственного числа А1 к минимальному ненулевому (положительному) А^. Так как процедура регуляризации матрицы [КО] фактически превращает три нулевых собственных числа матрицы [КО] в ненулевые (причём их можно изменить в ту или другую сторону), и так

как изменяя их, можно изменить число обусловленности с([К]), то нетрудно понять, что минимальным число с([К]) будет тогда, когда эти три собственных числа попадут внутрь спектра, то есть не будут ни максимальным, ни минимальным собственными числами.

Таким образом, разумно оценить обусловленность матрицы [КО] числом с0([КО]), а для его вычисления или оценки, в свою очередь, нужно оценить максимальное и минимальное положительные собственные числа матрицы [КО].

Так как матрица [КО] является суммой по всем КЭ неотрицательно-определённых матриц, каждая из которых построена по локальной матрице жёсткости некоторого КЭ, то было проведено исследование зависимости максимального и минимального собственных чисел матрицы такой частичной суммы в зависимости от различного выбора слагаемых.

Результаты вычислений для исходной сетки приведены на рис. 1, 2.

Область предполагалась имеющей прямоугольную форму. Исходная элементная сетка была образована взаимно перпендикулярными семействами линий, параллельных координатным осям, причём линии одного семейства располагались на равном расстоянии друг от друга, а каждый полученный прямоугольник был разделён однотипной диагональю на два треугольных КЭ. Всего имелось 20 прямоугольников в горизонтальном ряду и 4 в вертикальном. Первоначально брался угловой прямоугольник и к нему последовательно в направлении оси ОХ присоединялись соседние прямоугольники до завершающего этот горизонтальный ряд ещё одного углового прямоугольника включительно (кривая 1 на рис. 1 для А1 и кривая 1 на рис. 2 для АД Затем к полученному ряду добавлялся сразу весь следующий горизонтальный ряд (точка 1 на рис. 1 для А1), затем ещё следующий (точка 2 на рис. 1 для А1) и последний (самая правая точка кривой 2 на рис. 1). Потом вычисления были проделаны аналогич-

ным образом в направлении оси ОУ (кривая 2 на рис. 1 для А1 и кривая 2 на рис. 2 для А/). Отметим, что точки для А/,, соот-

ветствующие точкам 1 и 2 на рис. 1, совпадают с самой правой точкой кривых 1 и 2 на рис. 2.

5,00Е+012 и 4.50Е+012 4,00Е+012 3,50Е+012 -| А| 3,00Е+012 2,50Е+012 ■ 2,00Е+012-1.50Е+012 1,00Е+012

■ Кривая 1 • Точка 1 АТочка 2 ▼ Кривая 2

10

N

15

20

Рис. 1. Максимальное собственное число А1. По оси ОХ отложено количество N прямоугольных элементов по горизонтали в частичной сумме. Кривая 1 соответствует прямоугольникам N^1, точка 1 — прямоугольнику 20*2, точка 2 — прямоугольнику 20*3, кривая 2 — прямоугольникам N^4, то есть суммированию вертикальных рядов

3.50Е+011 3.00Е+011 2.50Е+011 2.00Е+011 -1.50Е+011 -1.00Е+011 -5.00Е+010-О.ООЕ+ООО

▼ ▼ ▼ ▼

...

-1—

10

Л"

15

■ Кривая 1 ▼ Кривая 2

—I—

20

Рис. 2. Минимальное ненулевое собственное число А/. По оси ОХ отложено количество N прямоугольных элементов по горизонтали в частичной сумме. Кривая 1 соответствует прямоугольникам ^1, кривая 2 — прямоугольникам ^4. Минимальное число обусловленности при 20*4 прямоугольниках имеет величину примерно 2,510

Анализ полученных результатов и основные выводы

Из проведённых вычислений следует:

1) при последовательном добавлении прямоугольников одного ряда максимальное собственное число монотонно возрастает, причём скорость его возрастания уменьшается, а минимальное собственное число убывает, причём скорость его убывания тоже уменьшается;

2) при добавлении следующего ряда, если такое добавление приводит к приближению отношения сторон полученного большого прямоугольника к 1:1 (то есть к приближению полученного большого прямоугольника к квадрату), то минимальное собственное число возрастает, если же отношение сторон отходит от 1:1, то минимальное собственное число уменьшается. Максимальное собственное число всегда увеличивается (что объясняется теоретически в теории матриц);

3) если первый ряд построить по более длинной стороне исходной области, а затем добавлять параллельные ему ряды КЭ, то минимум минимального числа обусловленности будет достигнут при завершении формирования первого ряда, а максимум максимального собственного числа — при завершении построения матрицы [КО].

Искажение сетки в процессе деформации приводит (непосредственно перед образованием первой трещины) к уменьшению минимума минимального собственного числа до 7,4-10 и увеличению максимума максимального собственного числа до 4,6-1012

Так как наличие сравнительно малых собственных чисел «разбалтывает» (выражение из [2]) решение в направлении собственных векторов, соответствующих малым собственным числам, так как такое «разбалтывание» в направлении указанных векторов нулевого собственного числа не влияет на НДС заготовки, и так как минимальное положительное собственное число достаточно велико даже при максимально искажённых КЭ, то в данном случае (для сетки 4 на 20 прямо-

угольников, разделённых каждый на два треугольника) можно считать задачу достаточно хорошо обусловленной, а результаты расчёта — достоверными.

Необходимо заметить, что при уменьшении КЭ с сохранением формы в п раз выражения Ъи Ъ., Ът, с., сс ст также уменьшатся в п раз, а площадь КЭ уменьшится в п2 раз. Следовательно, матрица [5] градиентов перемещений увеличится в п раз, а локальная матрица жёсткости не изменится. Из этого следует, что при использовании более мелких КЭ из-за увеличения количества КЭ произойдёт увеличение максимального и уменьшение минимального ненулевого собственного числа глобальной матрицы жёсткости. Таким образом, при измельчении сетки рано или поздно минимальное собственное число станет достаточно малым для того, чтобы «разбалтывание» решения СЛАУ в направлении соответствующего собственного вектора стало неприемлемым. В этот момент задача решения данной СЛАУ станет некорректной и для получения вектора придётся применять методику из [1].

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

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

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

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

Библиографический список

1. Тихонов, А.Н. Методы решения некорректных задач [Текст]/ А.Н. Тихо-

нов, В.Я. Арсенин. — М.: Наука, 1986. — 288 с.

2. Фаддеев, Д.К. Вычислительные методы линейной алгебры [Текст]/, Д.К. Фаддеев, В.Н. Фаддеева. — СПб.: Лань. 2009. — 736 с.

3. Chang, T. M., Swift H. W. Shearing of metal bars [Текст] / T. M. Chang, H. W. Swift // Journal of the Institute of Metals // vol. 78, 1950. - 119-146p.

NUMERICAL SIMULATION OF SEPARATION PROCESSES IN MECHANICAL WORKING

© 2012 K. M. Ivanov1, P. M. Vinnik1, V. N. Ivanov2

1Baltic State Technical University «Voenmekh» named after D. F. Ustinov 2Saint-Petersburg Institute of Machine Building

The paper deals with the methodological peculiarities of numerical simulation of separation processes in the context of the finite element method. It is shown that a poorly conditioned global stiffness matrix reflects the essence of the separation process. A method for monitoring the reliability and accuracy of the calculated data by means of the spectrum of a non-regularized global stiffness matrix is put forward. The possibility of assessing the spectrum of a global stiffness matrix through the sums of local matrices is shown.

Finite element method, separation processes, stiffness matrices, conditionality number.

Информация об авторах

Иванов Константин Михайлович, доктор технических наук, профессор, ректор, заведующий кафедрой «Технология производства артиллерийских систем и боеприпасов», Балтийский государственный технический университет «Военмех» им. Д.Ф. Устинова. E-mail: [email protected]. Область научных интересов: теория технологии обработки давлением, нелинейная механика и общая технология.

Винник Петр Михайлович, кандидат физико-математических наук, доцент кафедры высшей математики, Балтийский государственный технический университет «Военмех» им. Д.Ф. Устинова. E-mail: [email protected]. Область научных интересов: разделительные процессы обработки металлов давлением, метод конечных элементов, матричные вычисления.

Иванов Владимир Николаевич, кандидат технических наук, доцент, профессор кафедры «Машины и технологии обработки металлов давлением», Санкт-Петербургский институт машиностроения (ЛМЗ-ВТУЗ). E-mail:

[email protected]. Область научных интересов: порошковая металлургия, нано-технологии.

Ivanov Konstantin Mikhailovich, doctor of engineering, professor, rector, head of department, Baltic State Technical University «Voenmekh» named after D. F. Ustinov. E-mail: [email protected]. Area of research: theory of mechanical working technology, nonlinear mechanics and general production engineering.

Vinnik Pyotr Mikhailovich, candidate of physics and mathematics, associate professor of the department of higher mathematics, Baltic State Technical University «Voenmekh» named after D. F. Ustinov. E-mail: [email protected]. Area of research: separation processes in plastic working of metals, finite element method, matrix computations.

Ivanov Vladimir Nikolaevich, candidate of engineering, professor of the department "Machines and technologies of plastic working of metals", Saint-Petersburg Institute of Machine Building. E-mail: [email protected]. Area of research: powder metallurgy, nanotechnology.