сотрудник
Брянск, Брянская область, Россия
сотрудник
Брянск, Брянская область, Россия
сотрудник
Брянск, Брянская область, Россия
УДК 519.63 Численные методы решения дифференциальных уравнений с частными производными
ББК 3297 Вычислительная техника
Работа посвящена получению матрицы жесткости высокоточного, плоского конечного элемента с 6 степенями своды в узле, для решения плоских задач теории упругости методом конечных элементов. В научной литературе представлены подобные высокоточные элементы высоких порядков. Однако, теоретические результаты этих работ достаточно далеки от их практического применения. В настоящей работе приводиться исчерпывающе подробный вывод матрицы жесткости высокоточного конечного элемента. По аналогии, может быть получена матрица жесткости тетраэдрального конечного элемента с 12 степенями свободы в узле. Для тестирования полученной матрицы жесткости, была написана программа, основанная на MKЭ, с помощью которой выполнен расчет консольной балки. Погрешность расчета перемещений составила всего 0,22 %. Вывод: представленная в статье матрица жесткости, с большим успехом может использоваться в численных методах расчета напряженно-деформированного состояния.
задача, теория, упругость, перемещения, деформация, напряжение, расчет, консольная балка
Введение
В настоящее время в методе конечных элементов (МКЭ) используется большое количество самых разнообразных конечных элементов. Наиболее перспективными для решения плоских задач теории упругости следует считать треугольные конечные элементы [2]. Во-первых, эти элементы позволяют более гибко производить конечно-элементную дискретизацию исследуемой области. Во-вторых, треугольная область обладает определенными преимуществами с точки зрения математической задачи двухмерной интерполяции [3].
На рис. 1 показан высокоточный конечный элемент с 6 степенями свободы в узле, в котором, в качестве степеней свободы используются перемещения и производные от перемещений. Такой конечный элемент будет обладать повышенной точность, поскольку связывает условиями непрерывности не только поля перемещений, но и поля деформаций.
Принятие дифференциала поля перемещений в качестве степени свободы упрощает расчет напряжений в узлах, поскольку компоненты тензора напряжений в узле выражаются через первые производные поля перемещений. По этой же причине имеется возможность задавать граничные условия в напряжениях. Настоящая статья посвящена выводу матрицы жесткости треугольного конечного элемента с шестью степенями свободы в узле.
|
Рис. 1. Треугольный конечный элемент, с шестью степенями свободы в узле Fig. 1. A triangular finite element with six degrees of freedom at the node |
Поля перемещений
Обозначим узлы треугольного конечного элемента (рис 2) буквами 
, 
и 
. Каждому узлу треугольника соответствуют координаты 
, 
и 
.
Для треугольного конечного элемента более естественно использование 
-координат [4]. В этом случае, поле перемещений внутри конечного элемента можно описать с помощью пары однородных кубических полиномов:
(1)
или более коротко:
. (2)
Рис. 2. Система координат 
Fig. 2. Coordinate system
Для того, чтобы в качестве степеней свободы использовать производные от перемещений вдоль сторон 
треугольника, выполним дифференцирование зависимостей (1).
. (3)
Видим, что дифференцирование векторной функции 
свелось к дифференцированию векторной функции 
. Найдем соответствующие производные. Для чего, каждой стороне треугольника 
поставим в соответствие координату 
с началом в младшем узле (рис. 2).
Рассмотрим сторону 
треугольника. Связь между координатой 
и декартовой системой 
задается очевидными соотношениями:
. (4)
Для других сторон могут быть выписаны аналогичные соотношения. Тогда производная сложной функции по координате 
запишем следующим образом:
. (5)
Функции 
входящие в (5), являются сложными функциями от координат 
. Поэтому:
. (6)
Имея в виду известные зависимости (4):
.
Для треугольника известны соотношения
,
учитывая которые, для различных сочетаний 
и 
из набора , , :
.
Подставляя эти зависимости в (9), получаем производные по конкретным переменным 
, 
, и 
:
.
Эти зависимости можно представить как произведение некоторой матрицы 
, строение которой очевидно, на вектор 
:
, где 
. (7)
Где:
Подставив эти зависимости в (3) получаем:
, где . (8)
Векторы узловых перемещений
Введем вектор узловых перемещений [5]. С этой целью, выполним некоторые преобразования (8). Умножим (8) на длину стороны треугольника 
. Получим:
. (9)
Введем определение производной от функции вдоль некоторой стороны треугольника .
.
В этом случае, зависимость (9) примет вид:
. (10)
Новое понятие производной обеспечивает независимость компонент матрицы 
от размеров стороны треугольника и делает зависимость (10) универсальной для любого конечного элемента при дифференцировании вдоль любой из его сторон.
Теперь можно определить векторы узловых перемещений. Введем два локальных вектора узловых перемещений с компонентами, упорядоченными вдоль направлений 
и 
.
. (11)
Индекс 
здесь указывает на то, что производные взяты вдоль соответствующих сторон треугольника по формуле (10). Вместе эти два вектора образуют глобальный вектор перемещений для треугольного конечного элемента, с компонентами, упорядоченными по направлениям и .
.
Введем также три локальных вектора узловых перемещений для узла треугольного конечного элемента
. (11)
где 
. Вместе эти три вектора образуют глобальный вектор узловых перемещений для конечного элемента, компоненты которого упорядочены по узлам
.
Между векторами 
и 
установим связь с помощью матрицы 
:
|
25 |
. (12)
Строение матрицы очевидно.
Вектор узловых перемещений с производными вдоль сторон треугольного конечного элемента мы будем использовать для определения компонент матрицы . Использовать же этот вектор для вывода матрицы жесткости и расчете напряжений нельзя. В этом случае необходим вектор перемещений, компонентами которого будут являться производные по аргументам и . Введем следующее определение производных от функций вдоль координатных осей и .
. (15)
В соответствии с (15) определим три локальных вектора узловых перемещений для узла треугольного конечного элемента:
. (16)
где: 
. Индекс в (16) указывает на то, что производные взяты вдоль координатных осей и .
Из (9):
.
Теперь можно установить связь между векторами 
из (11) и 
из (16) с помощью матрицы 
:
(17)
где:
.
Здесь 
нулевая матрица размером 6×6.
Определим коэффициенты 
, 
(где 
) интерполяционных полиномов (1) Первые 9 коэффициентов определим из условий непрерывности перемещений и производных вдоль сторон треугольника в его узлах.
Для определения коэффициентов 
и 
необходим ещё один узел. Этот узел можно задать в центре тяжести треугольного элемента, что нежелательно. Существует много способов доопределения «лишних» неизвестных, которые довольно подробно излагаются в литературе. В результате, можно получить матрицу преобразования [Z], связывающую компоненты вектора перемещения внутри конечного элемента с узловыми перемещениями следующей зависимостью:
.
Где 
, 
, 
.
Введем матрицу функций формы:
.
Тогда окончательно:
. (18)
С целью облегчения расчета тензора напряжений, перейдем к использованию вектора узловых перемещений 
, компонентами которого являются производные от перемещений вдоль координатных осей и , для этого, подставим зависимость (17) в (18):
. (19)
Деформации и напряжения
Введем вектор деформаций:
.
и запишем уравнения Коши в виде:
, (20)
Подставим (19) в (20):
.
В результате действия матрицы 
на матрицу 
будет порождена матрица градиентов 
:
. (21)
Уравнение (21) дает искомую зависимость между деформациями и узловыми перемещениями.
Введем вектор напряжений:
.
и запишем связь между напряжениями и деформациями в виде:
, (22)
где 
- матрица упругости размером 
компоненты, которой определяются законом Гука Искомую зависимость получим, если в (22) подставить уравнения (21):
. (23)
Система разрешающих уравнений МКЭ
Введем обозначения:
- работа внешних сил
- потенциальная энергия деформированного тела
- энергия системы внешних и внутренних сил
Для вывода разрешающих уравнений МКЭ воспользуемся формулой Лагранжа [6]:
.
. (24)
Потенциальная энергия определяется формулой Клайперона [6]:
. (25)
Найдем потенциальную энергию 
для ко нечного элемента с номером 
, подставив в (25) уравнения (21) и (23):
. (26)
Работа внешних сил для узла конечного элемента с номером :
,
где: 
.
С помощью матрицы 
, введем в зависимость (26) вектор из (16)
,
где:
.
|
27 |
:
.
Выполнив суммирование по всем конечным элементам, получим:
.
Минимизируя полученный функционал приходим к разрешающей системе уравнений МКЭ:
. (27)
Матрица жесткости конечного элемента
Интеграл в левой части (27) представляет собой матрицу жесткости конечного элемента 
. Считая, что толщина конечного элемента изменяется линейно, запишем:
(L)=hiLi +hjLj+hkLk.
где: 
- толщина конечного элемента в узлах 
и учитывая, что компоненты матрицы 
являются постоянными величинами, получаем матрицу жесткости конечного элемента:
. (28)
.
Тестовый расчет
Используя полученные зависимости, была написана программа для тестирования матрицы жесткости конечного элемента и выполнен расчет перемещений для консольной балки (Рис. 3). Расчет имеет чисто математический характер, поэтому величины и размерности физических величин не приводятся.
Рис. 3. Триангуляция консольной балки для выполнения тестового расчета
Fig. 3. Triangulation of a cantilever beam for performing a test calculation
Расчет перемещений для точки балки (Рис. 3), по точным формулам, составил величину 2.744, а при расчете по методу конечных элементов получено перемещение 2.738. Таким образом, погрешность расчета составила величину 0.22%.
Заключение
В представленной работе подробно представлен вывод матрицы жесткости высокоточного конечного элемента с 6 степенями свободы в узле, для решения плоских задач теории упругости. Тестовый расчет подтвердил высокую точность конечного элемента. По полной аналогии, для решения объемных задач, можно получить матрицу жесткости тетраэдрального конечного элемента с 12 степенями свободу в узле.
1. Корнеев В.Г. Схемы метода конечных элементов высоких порядков точности. – Л.: Изд-во Ленингр. ун-та, 1977. 208 с.
2. Галлагер Р. Метод конечных элементов. Основы. – М.: Мир, 1984. – 428 с.
3. Сегерлинд Л. Применение метода конечных элементов. – М.: Мир, 1979. – 392 с.
4. Zienkiewicz O.C., Taylor R.L. The Finite Element Method. Vol. 1: The Basis. – Butterworth Heinemann, 2000. – 707p.
5. Баландин, М.Ю., Шурина Э.П. Векторный метод конечных элементов: Учеб. пособие. — Новосибирск: Изд-во НГТУ, 2001. — 69 с.
6. Ern A., Guermond J.L. Theory and practice of finite elements. Applied Mathematical Sciences. 2004;159.
7. Zienkiewicz O.C., Taylor R.L. The Finite Element Method. Vol. 2: Solid Mechanics. – Butterworth Heinemann, 2000. – 459 p.
8. Shames I.H., Cozzareli F.A. Elastic and Inelastic Stress Analysis, revised edition. – Washington: Taylor & Francis, DC, 1997. – 187 p.
9. Di P.D.A., Ern A. Mathematical aspects of discontinuous Galerkin methods. Mathématiques et Applications. 2012;69.
10. Стринг Г., Фикс Дж. Теория метода конечных элементов. – М.: Мир, 1977. – 350 с.
