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

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

3.1 Скорости перемещения

Пусть скорости перемещения в любой точке внутри элемента задаются вектор-столбцом:

(22)

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

(23)

(24)

(25)

где - скорость точки в плоскости заданной точками ;

- координата точки в плоскости заданной точками ;

-функции положения определяемые так же как и в 2D случае, а именно уравнениями (9-11);

- скорости по оси элементов соответственно;

- координаты узлов - соответственно (рис. 5).

Рис. 5. Точка O(x,y,z) в призматическом конечном элементе

В данном случае и так как данные функции скоростей не зависят от .

Уравнение (22) в более полной матричной форме запишется:

(26)

3.2 Скорости деформации

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

Пользуясь уравнением (22), скорости деформаци в матричной форме могут быть записаны в виде:

математический моделирование формоизменение сжатие

(27)

(28)

Матрица тензора малых деформаций:

(29)

Интенсивность скорости деформации:

(30)

Значения компонентов матрицы тензора малых деформаций ( из уравнения (24), из уравненя (25)):

(31.1)

(31.2)

(31.3)

(31.4)

(31.5)

(31.6)

3.3 Матрица жёсткости

Матрица из уравнения (28) состоит из трёх подматриц , каждая из которых имеет размерность . Зная вид функций формы из уравнения (9), можно определить компоненты каждой из подматриц:

(32)

(33)

(34)

Связь между напряжениями и скоростями деформации , характеризующая физическое состояние среды описывается соотношением:

(35)

где - накопившиеся напряжения в материале к рассматриваемому моменту времени, определяемые с помощью анализа истории нагружения материала. Уравнение (35) описывается подробнее таким образом:

(36)

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

(37)

Обзначим , тогда матрица и вектор , исходя из уравнений (35-36), имеют вид:

(38)

Матрица жёсткости -го элемента (размерность 9x9) имеет вид:

(39)

где каждый из компонентов - подматрица, имеющая размерность 3x3:

(40)

На первом этапе вычисляются компоненты матрицы - это проинтегрированная по от до из уравнения (25):

(41.1)

(41.2)

(41.3)

(41.4)

(41.5)

(41.6)

(41.7)

(41.8)

(41.9)

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

На втором этапе производится интегрирование результатов полученных выше по поверхности треугольника ( из уравнения (25)). Для линейных функций резальтат интегрирования - значение функции в центре ( треугольника помноженное на его площадь :

(42.1)

(42.2)

(42.3)

(42.4)

(42.5)

(42.6)

(42.7)

(42.8)

- нелинейный компонент, следовательно, для интегрирования по треугольнику, необходимо проинтегрировать функцию по и по . В данном случае удобно перейти из декартовой системы координат к L-кородинатам, которые являются функцями положения , каждая из которых изменяется от 0 до 1, при чём , так как :

(43)

Теперь интегрируется сначала по от до , после чего по от 0 до 1 (рис. 6):

(44)

Рис. 6. Угол

Пример вычисления для :

(45)

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

(46)

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

(47)

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

(48)

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

Итого :

(49)

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

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

Рис. 7. Разбиенкие треугольника рекуррентной функцией

3.4 Вектор правой части

Узловые силы, входящие в правую часть: {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 мм по оси и бесконечно длинный по оси ).

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