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

При рассмотрении твёрдого тела, подвергающегося механической деформации, задача математического моделирования формоизменения твёрдого тела заключается в построении вектор функции описывающей скорости частиц внутри его объёма [14]. В рассматриваемой задаче на границе твёрдого тела могут быть заданы следующие три типа граничных условий (рис. 2):

Рис. 2. Граничные условия

1)

2)

3)

(1)

где - компоненты тензора напряжений;

- компоненты нормали к поверхности в рассматриваемой точке;

- фрагменты контура тела, на которых заданы граничные условия,

- вектор распределённых поверхностных сил, действующих на фрагмент границы

, - единичные вектора, соответствующие ортам декартовой системы координат,

- компоненты скорости перемещения частиц среды, - скорости перемещения, заданные на фрагмент границы .

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

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

При такой постановке решается квазистатическая задача, то есть моделирование осуществляется итерационно, а именно: весь временной интервал, в котором решается задача, разбивается на малые подыинтервалы времени ?t. На каждом из подынтервалов считается, что скорость перемещений не изменяется, соответственно на каждом шаге выполняется уравнение равновесия:

(2)

где - напряжение.

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

(3)

где - компоненты тензора напряжений;

- скорости перемещения частиц среды;

- кажущийся коэффициент вязкости;

- символ Кронекера:

(4)

- накопленное гидростатическое напряжение:

(5)

где - коэффициент объёмного сжатия;

- относительное изменение объёма.

Итого, система уравнений, описывающая формоизменение материала:

6)

Глава 2. Конечно элементное решение и проблемы его классического применения

Узловые силы, входящие в правую часть: {R} - сосредоточенные внешне силы (в нашем случае = 0); - действие равномерно распределённых нагрузок интенсивности (в нашем случае = 0); - действие массовых сил (в нашем случае = 0); - вклад гидростатического давления:

(50)

где введены следующие обозначения:

Компоненты вектора правой части до интегрирования не зависят от Z, значит после интегрирования по Z, от 0 до из уравнения (25), будет:

(51)

Компоненты вектора правой части после интегрирования по Z зависят от X, Y линейно, значит результат их интегрирования по треугольнику : значение функции, проинтегрированной по Z, в центре треугольника (, помноженное на его площадь :

(52)

Глава 4. Результаты

4.1 Проверка полученных формул в сравнении с 2D случаем

Призматический элемент с координатами соответствует в геометрическом смысле треугольному элементу в 2D случае (с высотой ), когда координаты обоих элементов идентичны. Значения компонентов в полученной матрице жёсткости для описанного выше призматического элемента (уравнение (40)), должны быть идентичны значениям компонентов полученным в матрице жёсткости для описанного выше треугольного элемента в 2D случае (уравнение(20)).

В рамках работы были написаны следующие программы:

1) программа, реализующая вычисления с помощью метода конечных элементов в 2D случае для треугольного элемента (описаны в главе 2);

2) программа, реализующая вычисления для призматического элемента (описаны в главе 3).

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

Таблица 1

Координаты узлов рассматриваемых элементов и свойства среды

Демонстрация работы программы в консольном окне для вычисления матрицы жёсткости треугольного элемента:

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

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

Таблица 2

Матрица жёсткости для треугольного элемента

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

Демонстрация работы программы в консольном окне для вычисления матрицы жёсткости призматического элемента (причём разбиение для расчёта компонента происходит на (рис. 7)):

В результате создан файл «K = 20; deltaT = 1; mu = 0.03.cvc» в который записана следующая матрица жёсткости:

Таблица 3

Матрица жёсткости для призматического элемента

Компоненты в обеих матрицах жёсткости совпадают (в таблице 3 зелёным цветом выделены компоненты, соответствующие компонентам с которыми они совпадают из таблицы 2; оранжевым цветом выделены компоненты вычисленные при разбиение на глубину 3).

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

Изучено влияние глубины разбиения (рис. 7) на время вычисления и значение компонента матрицы жёсткости на примере призматического элемента, где - разные:

Таблица 4

Координаты узлов призматического элемента и свойства среды

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

Рис. 10 График зависимости значений компонента от глубины разбиения

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

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

Рис. 11. График зависимости времени вычисления компонента от глубины разбиения

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

4.3 Реализация тестовых расчетов по осадке образца кубической формы

Алгоритмы, описанные в главе 3, были реализованы в форме программных компонентов системы имитационного моделирования процессов обработки материалов давлением.

Проверка работы программы проводилась на примере осадки образца кубической формы (размеры которого: 20х20х20 мм) с помощью пресса (размеры которого: 20 мм по оси и бесконечно длинный по оси ).

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

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

Рис. 12. Триангуляция фрагмента образца

Испытание проводится до тех пор, пока правая стенка образца не достигнет значения 12 мм по оси X.

Рис. 13. Ввод параметров процесса

Значения по оси Z испытуемого образца в результате испытания:

Рис. 14. Результат испытания: значения по оси Z по осадке образца кубической формы

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

Сужение исследуемого фрагмента образца по оси Y произошло с до . Объём образца в начале испытания:

(53)

Объём образца в конце испытания:

(54)

Происходит допустимая потеря объёма, из-за накопленных ошибок компьютерного округления, а именно:

(55)

Следовательно, выполняется условие постоянства объёма при пластической деформации: [16].

Значение накопленной деформации [17]:

(56)

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

Рис. 15. Результат испытания: накопленная деформация

4.4 Реализация расчетов по моделированию испытания на сжатие с плоской деформацией

Программные компоненты системы имитационного моделирования были применены для исследования испытания на сжатие с плоской деформацией. Испытания проводились на образце, размеры которого: 15 мм по оси , 10 мм по оси , 20 мм по оси . Размеры пресса: 5 мм по оси , 20 мм по оси .

Сначала проводится триангуляция, то есть генерируется сетка образца:

Рис. 16. Триангуляция фрагмента образца

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

Рис. 17. Ввод параметров процесса

В результате испытания, сетка выглядит следующим образом:

Рис. 18. Результат испытания: сетка

Значения по оси Z испытуемого образца:

Рис. 19. Результат испытания: значения по оси Z

Распределение скоростей деформации в образце:

Рис. 20. Результат испытания: скорости деформации

Значения накопленной деформации:

Рис. 21. Результат испытания: накопленная деформация

Заключение