Интерполяционные полиномы для дискретизированной области

При обсуждении в предыдущем разделе интерполяционных соотношений для отдельных конечных элементов (симплекс – элементов) мы не фиксировали числовые значения узловых координат и ориентацию КЭ, выбирая их так, как было удобно для изложения сути проблемы. Подобный «произвол» весьма является важным достоинством метода КЭ, поскольку свобода выбора размеров (при выборе узловых координат) и ориентации КЭ позволяет составлять самые общие вычислительные алгоритмы и подпрограммы, моделирующие поведение отдельных КЭ. Обычно эти подпрограммы составляют основу библиотек конечных элементов в САПР теплового, прочностного и других видов анализа конструкции. Указанные подпрограммы могут быть далее использованы без изменения при рассмотрении областей с самыми разнообразными границами. Поскольку при решении задачи анализа поведения ИТО в заданной области производится ее дискретизации (разбиение на КЭ), то будем называть указанную область – дискретизированной областью (D-область).

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

Систему координат, связанную с элементом, будем называть местной, а систему координат, в которой задана D-область – глобальной.

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

Рассмотрим механизм преобразования интерполяционных соотношений, записанных в глобальной системе координат, к виду, представляющему эти соотношения в глобальной системе координат. С этой целью рассмотрим треугольный элемент, в котором скалярная величина представляется согласно как: j = Ni Фi + Nj Фj + Nk Фk.

Поместим местную систему координат в центре элемента, имеющего координаты (Xc, Yc), где:

Обозначив через s (t) – абсциссу (ординату) местной системы координат, запишем формулы преобразования координат (рисунок 3.2.1):

x = Xc + s; y = Xc + t

Функция формы Ni в глобальной системе координат, как было установлено ранее, имеет вид:

Ni = 0, 5 А ^–1 [Ai + Bi x + Ci y ]

Рис. 3.2.1

Подставляя сюда вместо x и y их выражения через s и t, получим:

Ni = 0, 5 А ^–1 [Ai + Bi (Xc+s) + Ci (Yc+t)] или:

Ni = 0, 5 А ^–1 [(Ai + Bi Xc + Ci Yc) + Bi s + Ci t ]

В результате преобразования Bi и Ci остаются неизменными и по-прежнему умножаются на независимые переменные. Константа же Ai – изменяется:

(Ai + Bi Xc + Ci Yc) = 2А/3

Таким образом, функция формы в местной системе координат равна:

Ni = 0, 5 А ^–1 [(2А/3) + (Yj –Yk) s + (Xk – Xj) t ]

Аналогично получим выражения для других функций формы:

Nj = 0, 5 А ^–1 [(2А/3) + (Yi –Yk) s + (Xk – Xi) t ]

Nk = 0, 5 А ^–1 [(2А/3) + (Yi –Yj) s + (Xj – Xi) t ]

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

где: R и R ^ – старая и новая области интегрирования, ê J÷ - отношение площадей в двух системах (J=Аxy/Ast). Форму элементов при переходе от местной системы координат к глобальной оставляют без изменений, поэтому R = R ^, кроме того, местная система и глобальная система координат являются декартовыми, поэтому J=1. Следовательно, в нашем случае:

Рассмотрим теперь задачу включения каждого конечного элемента, заданного в местной системе координат, в рассматриваемую D – область. Для этого необходимо выразить интерполяционные уравнения для каждого КЭ, используемого в ансамбле, через глобальные координаты и глобальные узловые значения.

С этой целью рассмотрим показанную на рисунке 3.2.2. пятиэлементную конфигурацию конечных двумерных симплекс – элементов, покрывающую некоторую D – область. Координаты всех узлов считаются известными. Узлы перенумерованы от 1 до 6, а порядковый номер конечного элемента указан на рисунке в скобках внутри соответствующего конечного элемента.

Условия Ф_b, выполняющиеся в узлах (b = 1, 2, …, 6), представляют собой глобальные степени свободы. Для составления интерполяционных полиномов, действующих внутри каждого конечного элемента, отметим символом «звездочка» один из узлов внутри каждого конечного элемента. Тогда, по аналогии с полученным ранее выражением, для элемента [1] можно записать:

j^[1] = N₂^[1]Ф₂ + N₃^[1]Ф₃ + N₁^[1]Ф₁

Здесь выражения N_q^[1] представляют функции формы конечного элемента [1] в q–м узле (q=1, 2, 3). Для их правильного вычисления следует подставлять значения глобальных координат. Только в этом случае выражение действительно позволит учесть соответствие индексов элемента [i, j, k] глобальным номерам узлов. Согласно принятой на рисунке нумерации узлов Э_q, соответствие [i® j® k] примет вид:

Э₁: [2®3®1], Э₂: [3®2®4], Э₃: [5®3®4], Э₄: [6®3®5], Э₅: [1®3®6]

Учитывая принятую нумерацию индексов, приходим к следующей совокупности уравнений для всех конечных элементов D – области:

Рис. 3.2.2

j^[1] = N₂^[1]Ф₂ + N₃^[1]Ф₃ + N₁^[1]Ф₁

j^[2] = N₃^[2]Ф₃ + N₂^[2]Ф₂ + N₄^[2]Ф₄

j^[3] = N₅^[3]Ф₅ + N₃^[3]Ф₃ + N₄^[3]Ф₄

j^[4] = N₆^[4]Ф₆ + N₃^[4]Ф₃ + N₅^[4]Ф₅

j^[5] = N₁^[5]Ф₁ + N₃^[5]Ф₃ + N₆^[5]Ф₆

Пример 3.2.1.

Получить функции формы N₆^[4] и N₆^[5] в заданной на рисунке 3.2.2 D – области.

Решение.

1. В соответствии с выражениями (9.11) и (9.9) запишем общее выражение для функции формы в произвольной точке элемента 4:

2. Пользуясь выражением (9.25), запишем соответствие индексов для конечного элемента Э₄ имеем: [i=6®j=3®k=5]. Заменяя индексы их значениями, получим для N₆^[4] :

3. Аналогичное соответствие индексов для конечного элемента Э₅ имеем: [i=1®j=3®k=6]. Заменяя индексы их значениями, получим для N₆^[5]:

Последние две формулы показывают, что функции формы N₆^[5] и N₆^[4] – совершенно разные функции, аппроксимирующие заданный функционал соответственно в пределах конечных элементов 5 и 4. Однако, в самом шестом узле обе эти функции принимают единичные значения, поскольку числители обоих формул в этой точке принимают значения определителя, равные 2А.