Введение В настоящее время в методе конечных элементов (МКЭ) используется большое количество самых разнообразных конечных элементов. Наиболее перспективными для решения плоских задач теории упругости следует считать треугольные конечные элементы [2]. Во-первых, эти элементы позволяют более гибко производить конечно-элементную дискретизацию исследуемой области. Во-вторых, треугольная область обладает определенными преимуществами с точки зрения математической задачи двухмерной интерполяции [3]. На рис. 1 показан высокоточный конечный элемент с 6 степенями свободы в узле, в котором, в качестве степеней свободы используются перемещения и производные от перемещений. Такой конечный элемент будет обладать повышенной точность, поскольку связывает условиями непрерывности не только поля перемещений, но и поля деформаций.Принятие дифференциала поля перемещений в качестве степени свободы упрощает расчет напряжений в узлах, поскольку компоненты тензора напряжений в узле выражаются через первые производные поля перемещений. По этой же причине имеется возможность задавать граничные условия в напряжениях. Настоящая статья посвящена выводу матрицы жесткости треугольного конечного элемента с шестью степенями свободы в узле.   Рис. 1. Треугольный конечный элемент, с шестью степенями свободы в узлеFig. 1. A triangular finite elementwith six degrees of freedom at the node  Поля перемещений Обозначим узлы треугольного конечного элемента (рис 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 : 25Us=E1Ups .                                                  (12)  Строение матрицы E1  очевидно.Вектор узловых перемещений с производными вдоль сторон треугольного конечного элемента мы будем использовать для определения компонент матрицы α10 . Использовать же этот вектор для вывода матрицы жесткости и расчете напряжений нельзя. В этом случае необходим вектор перемещений, компонентами которого будут являться производные по аргументам x  и y . Введем следующее определение производных от функций uL  вдоль координатных осей x  и y . uL,1=∂uL∂x, uL,2=∂uL∂y .                                        (15)В соответствии с (15) определим три локальных вектора узловых перемещений для узла p  треугольного конечного элемента:UpxT=UpUp,1Up,2VpVp,1Vp,2 .                                   (16)где: p∈i, j, k . Индекс x  в (16) указывает на то, что производные взяты вдоль координатных осей x  и y . Из (9):u(L),,r=∂u(L)∂Srlr=lr∂u(L)∂x∂x∂Sr+∂u(L)∂y∂y∂Sr     =lr∂u(L)∂xcr-∂u(L)∂ybr .Теперь можно установить связь между векторами Ups  из (11) и Upx  из (16) с помощью матрицы Mpsx :Ups=MpsxUpx                                                   (17)где:MpSx=MiSx060606MjSx060606MkSx .Здесь 06  нулевая матрица размером 6×6.MiSx=1000cj-bj0ck-bk    MjSx=1000ck-bk0ci-bi    MkSx=1000ci-bi0cj-bj   Определим коэффициенты αl , βl  (где l∈1, 2, …, 10 ) интерполяционных полиномов (1) Первые 9 коэффициентов определим из условий непрерывности перемещений и производных вдоль сторон треугольника в его узлах.Для определения коэффициентов α10  и β10  необходим ещё один узел. Этот узел можно задать в центре тяжести треугольного элемента, что нежелательно. Существует много способов доопределения «лишних» неизвестных, которые довольно подробно излагаются в литературе. В результате, можно получить матрицу преобразования [Z], связывающую компоненты вектора перемещения внутри конечного элемента с узловыми перемещениями следующей зависимостью: uL=ΩLZ00ZUs .Где        uL=&u(L)&v(L) ,   ΩL=L10T010T010TL10T ,   Us=&Us&Vs .Введем матрицу функций формы:NL=ΩLZ00Z .Тогда окончательно:uL=NLUps .                                                           (18)  С целью облегчения расчета тензора напряжений, перейдем к использованию вектора узловых перемещений Upx , компонентами которого являются производные от перемещений вдоль координатных осей x  и y , для этого, подставим зависимость (17) в (18): uL=NLMsxUpx .                                          (19)Деформации и напряженияВведем вектор деформаций:εT=εxτxyεy=ε11ε12ε22 .и запишем уравнения Коши в виде:εL=SuL ,                                                         (20)Подставим (19) в (20):εL=SNLMsxUpx .В результате действия матрицы S  на матрицу NL  будет порождена матрица градиентов BL : εL=BLMsxUpx .                                                  (21)Уравнение (21) дает искомую зависимость между деформациями и узловыми перемещениями. Введем вектор напряжений:σT=σxτxyσy=σ11σ12σ22 .и запишем связь между напряжениями и деформациями в виде:σL=DεL ,                                                     (22) где D  - матрица упругости размером 3×3  компоненты, которой определяются законом Гука Искомую зависимость получим, если в (22) подставить уравнения (21): σL=DBLMsxUpx .                                                     (23) Система разрешающих уравнений МКЭВведем обозначения:A  - работа внешних силΠ  - потенциальная энергия деформированного телаΛ  - энергия системы внешних и внутренних силДля вывода разрешающих уравнений МКЭ воспользуемся формулой Лагранжа [6]:δΛ=0 .Λ=Π-A .                                                         (24)Потенциальная энергия определяется формулой Клайперона [6]:Π=12VεTσdV .                                                  (25)Найдем потенциальную энергию Πν  для ко          нечного элемента с номером ν , подставив в (25) уравнения (21) и (23):Πν=12VUpxTMνsxTBνLTDνBνLMνsxUpxdV .                     (26)Работа внешних сил для узла p  конечного элемента с номером ν :Apν=UpνTPpν ,где: UpνT=UpνUpν; PpνT=PpxνPpyν .С помощью матрицы F2 , введем в зависимость (26) вектор Upx  из (16) Apν=UpxTF2Ppν ,где:F2=100000000100 . 27Тогда, для конечного элемента с номером ν : Aν=p=13UpνxTF2Ppν .Выполнив суммирование по всем конечным элементам, получим:Λ=12ν=1ΕVUpxTMνsxTBνLTDνBνLMνsxUpxdV-ν=1Εp=13UpνTF2Ppν .Минимизируя полученный функционал приходим к разрешающей системе уравнений МКЭ:ν=1ΕVMνsxTBνLTDνBνLMνsxUpxdV=ν=1Εp=13F2Ppν .                 (27) Матрица жесткости конечного элемента Интеграл в левой части (27) представляет собой матрицу жесткости конечного элемента k . Считая, что толщина конечного элемента изменяется линейно, запишем: (L)=hiLi +hjLj+hkLk. где: hi, hj, hk  - толщина конечного элемента в узлах i, j, k  и учитывая, что компоненты матрицы Mνsx  являются постоянными величинами, получаем матрицу жесткости конечного элемента:  kν=MνsxTVBνLTDνBνLdVMνsx .                                      (28).Тестовый расчет Используя полученные зависимости, была написана программа для тестирования матрицы жесткости конечного элемента и выполнен расчет перемещений для консольной балки (Рис. 3). Расчет имеет чисто математический характер, поэтому величины и размерности физических величин не приводятся.   Рис. 3. Триангуляция консольной балки для выполнения тестового расчетаFig. 3. Triangulation of a cantilever beam for performing a test calculation  Расчет перемещений для точки A  балки (Рис. 3), по точным формулам, составил величину 2.744, а при расчете по методу конечных элементов получено перемещение 2.738. Таким образом, погрешность расчета составила величину 0.22%.   Заключение В представленной работе подробно представлен вывод матрицы жесткости высокоточного конечного элемента с 6 степенями свободы в узле, для решения плоских задач теории упругости. Тестовый расчет подтвердил высокую точность конечного элемента. По полной аналогии, для решения объемных задач, можно получить матрицу жесткости тетраэдрального конечного элемента с 12 степенями свободу в узле.