Этап 3: объединение конечных элементов в ансамбль

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

Рассмотрим процедуру составления ансамбля конечных элементов для задачи нахождения поля температур в стержне (рис. 2.18, а). Кусочно-элементная модель области приведена на рис. 2.18, б, а функция отдельного элемента определяется уравнением (2.50).

Пример использования МКЭ для расчета одномерного температурного поля в однородном стержне длиной L и площадью поперечного сечения S (рис. 2.17). [1, 6]

Рис. 2.17. Однородный стержень, находящийся
под воздействием теплового потока

 

Один конец стержня закреплен, к нему подводится тепловой поток q заданной интенсивности. На свободном конце стержня происходит конвективный теплообмен с внешней средой. Коэффициент теплообмена α, температура окружающей среды Т_* известны. Вдоль боковой поверхности стержень теплоизолирован.

Температурное поле в стержне описывается уравнением теплопроводности, которое в одномерном приближении имеет вид

 

. (2.56)

 

Краевые условия определяются уравнениями:

 

при х= 0; (2.57, а)

при x = L. (2.57, б)

 

Искомое температурное поле является непрерывной функцией координаты х (рис. 2.18, а).

Рис. 2.18. Расчет одномерного температурного поля
в однородном стержне методом МКЭ.

 

В МКЭ стержень разбивается произвольным образом на конечные элементы (в данном случае отрезки неравной длины). На каждом элементе непрерывная функция Т (х) аппроксимируется некоторой линейной зависимостью, как показано на рис. 2.18, б (в скобках указаны номера элементов). Аппроксимирующая кусочно-линейная функция определяется через узловые значения T₁-Т₆, которые в общем случае сначала неизвестны и подлежат определению в МКЭ.

Запишем соответствие между произвольными номерами i, j, уравнения (2.50), и глобальными номерами узлов рассматриваемой дискретной модели для

элемента 1: i = 1, j = 2, (2.58, а)

элемента 2: i = 2, j = 3, (2.58, б)

элемента 3: i = 3, j = 4, (2.58, в)

элемента 4: i = 4, j = 5, (2.58, г)

элемента 5: i = 5, j = 6. (2.58, д)

 

Подставив значения номеров узлов (2.58) в (2.50), получим:

 

(2.59)

где верхние индексы в скобках относятся к номеру элемента.

В выражениях для функций формы элемента (2.49) значения произвольных номеров i, j также следует изменить в соответствии с (2.58). Тогда значения , определяются по формулам:

 

 

Очевидно, что и не равны друг другу даже в случае равенства длин элементов L ⁽²⁾ и L ⁽³⁾. При известных значениях узловых величин Т ₁- Т ₆ уравнения (2.59) позволяют определить значение температуры в любой точке стержня.

Рассмотрим пример объединения элементов двумерной области в ансамбль, который потребуется для иллюстрации дальнейших этапов МКЭ. Треугольная область разбита на элементы треугольной формы, как показано на рис. 2.19.

Для обозначения узлов отдельных элементов по-прежнему используются номера i, j, k, начиная с произвольного узла (на рисунке отмеченного звездочкой) против часовой стрелки. Соответствие между этими обозначениями и глобальными номерами узлов следующее:

 

элемента 1: i = 1, j = 2, k = 4, (2.60, а)

элемента 2: i = 1, j = 3, k = 5, (2.60, б)

элемента 3: i = 2, j = 5, k = 4, (2.60, в)

элемента 4: i = 4, j = 5, k = 6, (2.60, г)

 

Рис. 2.19. Пример составления ансамбля конечных элементов для двухмерной треугольной области

 

Подставляя значения (2.60) в (2.54), получим:

 

(2.61)

 

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

Система (2.61) — сокращенная форма математического описания модели. Расширенная форма имеет вид:

 

 

или в матричной форме

φ = NФ. (2.62)

В САПР с целью уменьшения объема памяти чаще используют сокращенную форму описания моделей (2.61). Расширенная форма описания моделей имеет некоторые преимущества при реализации следующих этапов алгоритма МКЭ.

Определение вектора узловых значений функций.

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

Пример минимизации функционала в задаче нахождении распределения температуры в однородном стержне (рис. 2.18). При указанном методе минимизируется функционал:

 

, (2.63)

где V — объем тела; S — площадь границы.

В функционал F входят оба граничных условия (2.57 а, б). При минимизации функционала используется множество функций элементов дискретизированной области. Для простоты вычислений будем считать, что стержень разбит всего на два элемента (в практических случаях этого недостаточно) (рис. 2.20).

 

Рис. 2.20. Пример минимизации функционала нахождения распределения температуры в стержне

 

 

Тогда

 

(2.64)

 

Функционал (2.63) удобно представить в виде

 

, (2.65)

где S₁ и S₂ — площади сечений стержня, на которых заданы граничные условия (2.57 а) и (2.57 б) соответственно.

 

Для вычисления объемного интеграла в (2.65) его необходимо разбить на два слагаемых в соответствии с принятой конечной элементной моделью:

 

. (2.66)

 

Производные в (2.66) вычисляются с учетом (2.64) и (2.49):

 

(2.67)

 

Подставив (2.67) в (2.66) и считая, что , получим

 

.

Второе и третье слагаемые в (2.65) вычисляются просто так как подынтегральным функциям соответствуют узловые значения Т ₁ и Т ₃:

 

;

,

где S ₁ и S ₂ — площади поверхностей, на которых заданы q и α (для рассматриваемого примера и ).

 

Значение функционала F вычисляется простым суммированием последних трех выражений:

 

.(2.68)

где и .

 

Для минимизации функционала F необходимо выполнение условий:

 

или в матричной форме

 

. (2.69)

В общем виде (2.69) можно представить так: , что соответствует (2.45).

Примечание. Матрица коэффициентов К в (2.69) по-прежнему называется матрицей жесткости, хотя по физическому смыслу данной задачи ее удобнее было бы назвать матрицей теплопроводности. Такое название матрицы К пришло из строительной механики, где МКЭ начал применяться раньше, чем в других областях техники.

Зная характеристики материала, из системы (2.69) можно определить узловые значения T ₁, Т ₂, Т ₃.

Из (2.68) и (2.69) нетрудно заметить, что однотипные конечные элементы вносят в эти выражения слагаемые одного вида. Поэтому при реализации МКЭ в САПР вклад элемента определенного типа в матрицу жесткости вычисляется только один раз, а затем используется во всех необходимых случаях. При этом алгоритм получения матрицы жесткости несколько отличается от описанного выше и состоит из следующих этапов:

Этап 1. Представление функционала F в виде суммы соответствующих функционалов для элементов.

Для рассмотренного примера , причем

 

;

.

Этап 2. Минимизация функционала каждого элемента отдельно (при этом вычисляются матрицы жесткости К ^(е) и векторы нагрузки В ^(е) для всех конечных элементов). В примере

.

Этап 3. Суммирование матриц жесткости и векторов нагрузки отдельных элементов [сумма приравнивается нулю, что позволяет получить систему КФ = В (2.45)].

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

Пример минимизации функционала в задаче о кручении стержня (рис. 2.21).

Рис. 2.21. Упругий стержень, находящийся под воздействием крутящего момента

 

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

 

.

Для определения напряжений необходимо знать форму сечения стержня. В случае квадратного сечения стержня из-за наличия симметрии при решении задачи достаточно рассмотреть 1/8 его часть.

Рассмотрение исследуемой части сечения стержня на элементы показано на рисунке 2.22.

 

Рис. 2.22. Поперечное сечение стержня, находящегося
под действием крутящего момента

 

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

 

.

 

Рассмотрим один из функционалов

 

 

и, введя обозначения

,

преобразуем выражение к виду

 

. (2.70)

 

С учетом (2.62) вектор g^(е) запишется в виде

 

, (2.71)

 

или g^(e) = D^(e)Ф, где D^(e) — матрица градиентов.

 

Подставив (2.71) в (2.70), получим

 

.

 

Минимизация функционала F ^(e) приводит к уравнению

 

,

где (2.72)

есть матрица жесткости, а

 

(2.73)

есть вектор нагрузки отдельного элемента.

Вычисление по (2.72) и (2.73) проведем на примере первого элемента. С учетом первого уравнения системы (2.61) матрицу градиентов D⁽¹⁾ можно представить в виде

 

где

 

Произведение матриц равно

,

где .

 

Все это выражение — постоянная величина, которая может быть вынесена из–под знака интеграла. С учетом сказанного и принимая для простоты толщину элемента равной 1, матрицу жесткости можно вычислить по формуле

 

. (2.74)

 

Для определения вектора нагрузки (2.73) удобно воспользоваться интегральной формулой, позволяющей вычислять интегралы по площади двумерных симплекс–элементов:

.

 

В частном случае при а = 1, b = 0, с = 0 формула имеет вид

 

. (2.75)

 

С учетом принятого допущения о единичной толщине элемента и формулы (2.75) вектор нагрузки В ⁽¹⁾ для первого элемента запишется в виде

 

.

 

Аналогичные вычисления необходимо проделать и для всех остальных элементов области. Принимая, что площади всех элементов равны, а свойства материала одинаковы, и объединяя уравнения для элементов, получим полную систему уравнений:

 

 

При реализации МКЭ в САПР форма (2.74) матрицы жесткости элемента неэффективна с точки зрения затрат ОП. Действительно, матрицы жесткости отдельных элементов имеют ту же размерность, что и глобальная матрица жесткости системы, а большинство элементов матрицы нулевые. В САПР с целью сокращения затрат ОП из матриц жесткости исключают нулевые элементы, строя их в сокращенной форме. Такой метод построения матриц называют методом прямой жесткости. При этом исключается необходимость хранения матриц большой размерности, но возникает потребность в специальной процедуре кодирования узлов элементов.

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

.

Вектор нагрузки также вычисляется в сокращенном виде:

.

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

Сочетание метода Галеркина с МКЭ приводит к системе уравнений [6]

 

где L(φ) — левая часть исходного дифференциального уравнения, описывающего непрерывную функцию φ. Технику получения разрешающей системы уравнений методом Галеркина легко проиллюстрировать на примере уже решенной выше задачи об отыскании температурного поля в однородном стержне (рис. 2.17), конечно–элементная модель которого представлена на рис. 2.21.

 

Применив метод Галеркина к (2.42), получим

 

. (2.76)

 

Подставим в (2.76) формулу дифференцирования произведений:

 

. (2.77)

 

Интерполяционная функция Т не сохраняет постоянства по длине стержня, поэтому интеграл в (2.77) можно представить суммой соответствующих интегралов для отдельных элементов. Так, второй интеграл в (2.77) можно представить в виде

 

.

 

Вычислим в этом уравнении интегралы, относящиеся к отдельным элементам:

; (2.78)

. (2.79)

 

С учетом (2.78) и (2.79)

….(2.80)

 

Первый интеграл в (2.77) на основании теоремы Остроградского–Гаусса преобразуется к виду

 

, (2.80)

где ; n — внешняя нормаль к рассматриваемой поверхности.

С учетом краевого условия (2.43 а) в точке для первого элемента интеграл (2.80) принимает вид

 

. (2.81)

 

Подставив значения , в (2.81), получим

 

. (2.82)

С учетом краевого условия (2.43 б) в точке второго элемента интеграл (2.82) запишется так:

…………..(2.83)

 

Просуммировав выражения вида (2.80) для первого и второго элементов и выражения (2.82) и (2.83) и приравняв сумму нулю, получим систему уравнений

 

, …...(2.84)

где .

 

Система (2.84) идентична системе (2.68), определенной путем минимизации соответствующим образом подобранного функционала.

Завершающим шагом этапа определения вектора силовых значений Ф является решение системы линейных алгебраических уравнений.

Примечание. Основные особенности этого шага — большая размерность и сильная разреженность матрицы коэффициентов системы. В связи с этим для реализации МКЭ в САПР разработаны специальные способы хранения матрицы жесткости, позволяющие уменьшить необходимый для этого объем ОП. Для нахождения узловых значений функций применяются методы преобразования и решения системы, направленные на снижение затрат машинного времени.

Нестационарные краевые задачи.

Во всех рассмотренных выше примерах МКЭ применялся для решения стационарных краевых задач. Алгоритм метода и особенности отдельных его этапов остаются неизменными при решении нестационарных задач, в уравнениях которых присутствуют не только частные производные по пространственным координатам, но и частные производные по времени.

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

 

. (2.85)

 

В этом случае при вариационной формулировке задачи минимизации подлежит функционал

. (2.86)

 

Функционал (2.86) отличается от минимизированного выше функционала (2.62) наличием дополнительного слагаемого в объемном интеграле. Вычислим вклад этого слагаемого в результирующую систему уравнений. Для этого обозначим

 

.

 

Принимая, что характеристики материала постоянны по всему объему стержня, перепишем данное уравнение в виде .

 

Интеграл F ₁ можно представить суммой интегралов, записанных для отдельных элементов:

. (2.87)

 

Используя соотношение T^(e) = N^(e)T^(e) и учитывая, что N ^(e) является функцией только координат, производную по времени в (2.87) представим в виде

 

.

 

Тогда для отдельного элемента функционал F ₁ записывается в форме

 

. (2.88)

 

Продифференцировав (2.88) по вектору Т ^(е) и приравняв результат нулю, получим

 

. (2.89)

 

Подставляя в объемный интеграл (2.89) соотношения для функций формы (2.50) и учитывая соотношение D^(e) = S^(e)dL^(e), найдем

(2.90)

 

После суммирования по элементам и с учетом (2.90) функционал F ₁ принимает вид

(2.91)

где H⁽¹⁾ = D⁽¹⁾S⁽¹⁾L⁽¹⁾/ 6 и H⁽²⁾ = D⁽²⁾S⁽²⁾L⁽²⁾/ 6.

 

Объединяя (2.91) с ранее полученным для той же задачи уравнением (2.68), результирующую систему уравнений для задачи (2.85) запишем в форме

. (2.92)

 

 

2.5.3. Программные комплексы конечно–элементного
анализа