Введение

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

На рис. 1 показан высокоточный конечный элемент с 6 степенями свободы в узле, в котором, в качестве степеней свободы используются перемещения и производные от перемещений. Такой конечный элемент будет обладать повышенной точность, поскольку связывает условиями непрерывности не только поля перемещений, но и поля деформаций.

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

Поля перемещений

Обозначим узлы треугольного конечного элемента (рис 2) буквами i , j  и k . Каждому узлу треугольника соответствуют координаты Xi, Yi , Xj, Yj  и Xk, Yk .

Для треугольного конечного элемента более естественно использование L  -координат [4]. В этом случае, поле перемещений внутри конечного элемента можно описать с помощью пары однородных кубических полиномов:

&uL=α1Li3+α2Lj3+α3Lk3+α4Li2Lj+α5Lj2Lk+α6Lk2Li+α7Lj2Li+α8Lj2Lk2+α9Li2Lk+α10LiLjLk&vL=β1Li3+β2Lj3+β3Lk3+β4Li2Lj+β5Lj2Lk+β6Lk2Li+β7Lj2Li+β8Lj2Lk2+β9Li2Lk+β10LiLjLk (1)

  или более коротко:

uL=α10⋅L10 .                                                                (2)

Рис. 2. Система координат Si, Sj, Sk

Fig. 2. Coordinate system Si, Sj, Sk

Для того, чтобы в качестве степеней свободы использовать производные от перемещений вдоль сторон Sr(r=i, j, k)  треугольника, выполним дифференцирование зависимостей (1).

∂uL∂Sr=∂∂Srα10L10=α10∂L10∂Sr .                                          (3)

Видим, что дифференцирование векторной функции uL  свелось к дифференцированию векторной функции L10 . Найдем соответствующие производные. Для чего, каждой стороне треугольника lrr=i, j, k  поставим в соответствие координату Sr  с началом в младшем узле (рис. 2).

Рассмотрим сторону jk  треугольника. Связь между координатой Si  и декартовой системой x, y  задается очевидными соотношениями:

&x=   ciliSi+Xj&y=-biliSi+Yj .                                                                  (4)

Для других сторон могут быть выписаны аналогичные соотношения. Тогда производная сложной функции по координате Sr  r=i, j, k  запишем следующим образом:

∂L10∂Sr=∂L10∂Li⋅∂Li∂Sr+∂L10∂Lj⋅∂Lj∂Sr+∂L10∂Lk ⋅∂Lk∂Sr .                                    (5)

Функции Lpp=i, j, k  входящие в (5), являются сложными функциями от координат x, y . Поэтому:

∂Lp∂Sr=∂Lp∂x⋅∂x∂Sr+∂Lp∂y⋅∂y∂Sr .                                                  (6)

Имея в виду известные зависимости (4):

∂Lp∂Sr=12⋅Δ⋅lrbpcr-brcp .

Для треугольника известны соотношения

cibk-ckbi=cjbi-cibj=ckbi-cjbk=ai+aj+ak=2Δ ,

учитывая которые, для различных сочетаний p  и r  из набора i , j , k :

∂Lp∂Sr=lr-1приr=i, p=k; r=j, p=i; r=k, p=j-lr-1приr=i, p=j; r=j, p=k; r=k, p=i0приr=p .

Подставляя эти зависимости в (9), получаем производные по конкретным переменным Si , Sj , и Sk :

&∂L10∂Si=1li-∂L10∂Lj+∂L10∂Lk&∂L10∂Sj=1lj-∂L10∂Lk+∂L10∂Li&∂L10∂Sk=1lk-∂L10∂Li+∂L10∂Lj .

Эти зависимости можно представить как произведение некоторой матрицы Trr=i, j, k , строение которой очевидно, на вектор L6 :

∂L10∂Sr=1lrTrL6 , где r∈i, j, k .                                              (7)

Где:

L6T=Li2, Lj2, Lk2, LiLj, LjLk, LkLi

Подставив эти зависимости в (3) получаем:

∂uL∂Sr=1lrα10TrL6 , где r∈i, j, k .                                         (8)

Векторы узловых перемещений

Введем вектор узловых перемещений [5]. С этой целью, выполним некоторые преобразования (8). Умножим (8) на длину стороны треугольника lr . Получим:

∂uL∂Sr⋅lr=α10TrL6 .                                                     (9)

Введем определение производной от функции uL  вдоль некоторой стороны треугольника lr .

ul,,r=∂uL∂Sr⋅lr .                                                             

В этом случае, зависимость (9) примет вид:

ul,,r=α10TrL6 .                                                  (10)

Новое понятие производной обеспечивает независимость компонент матрицы α10  от размеров стороны треугольника и делает зависимость (10) универсальной для любого конечного элемента при дифференцировании вдоль любой из его сторон.

Теперь можно определить векторы узловых перемещений. Введем два локальных вектора узловых перемещений с компонентами, упорядоченными вдоль направлений x  и y .

&UsT=UiUi,,jUi,,kUjUj,,kUj,,iUkUk,,iUk,,j&VsT=ViVi,,j Vi,,k Vj Vj,,k Vj,,iVk Vk,,i Vk,,j .                   (11)

Индекс s  здесь указывает на то, что производные взяты вдоль соответствующих сторон треугольника по формуле (10). Вместе эти два вектора образуют глобальный вектор перемещений для треугольного конечного элемента, с компонентами, упорядоченными по направлениям x  и y .

UsT=UsTVsT .

Введем также три локальных вектора узловых перемещений для узла p  треугольного конечного элемента

UpsT=UpUp,,rUp,,tVpVp,,rVp,,t .                                               (11)

где p∈i, j, k; r∈i, j, k; t∈i, j, k . Вместе эти три вектора образуют глобальный вектор узловых перемещений для конечного элемента, компоненты которого упорядочены по узлам

UpsT=UisTUjsTUksT .

Между векторами Us  и Ups  установим связь с помощью матрицы E1 :

Us=E1Ups .                                                  (12)

Строение матрицы E1  очевидно.

Вектор узловых перемещений с производными вдоль сторон треугольного конечного элемента мы будем использовать для определения компонент матрицы α10 . Использовать же этот вектор для вывода матрицы жесткости и расчете напряжений нельзя. В этом случае необходим вектор перемещений, компонентами которого будут являться производные по аргументам x  и y . Введем следующее определение производных от функций uL  вдоль координатных осей x  и y .

uL,1=∂uL∂x, uL,2=∂uL∂y .                                        (15)