Лекции - Модели и методы АПР
Подождите немного. Документ загружается.
Рис. 2. Расчет одномерного температурного поля в однородном стержне методом МКЭ.
В МКЭ стержень разбивается произвольным образом на конечные элементы, которые в
данном случае являются отрезками неравной длины. На каждом элементе непрерывная
функция Т(х) аппроксимируется некоторой линейной зависимостью, как показано на рис. 2, б (в
скобках указаны номера элементов). Аппроксимирующая кусочно-линейная функция
определяется через узловые значения T
1
-Т
6
, которые в общем случае сначала неизвестны и
подлежат определению в МКЭ.
В общем случае алгоритм МКЭ состоит из четырех этапов:
Этап 1. Выделение конечных элементов (разбиение заданной области на
конечные элементы).
Этап 2. Определение аппроксимирующей функции для каждого элемента
(определение функции элемента). На данном этапе значение непрерывной
функции
е
в произвольной точке е-го конечного элемента
аппроксимируется полиномом
0
ARA
eе
, где А
(е)
- вектор-строка
коэффициентов полинома; А
о
- свободный член; R=(х, у, z) - вектор координат
в рассматриваемой точке.
Задача этапа заключается в определении неизвестного вектора А
(е)
и
свободного члена А
о
. Для этого, используя условие непрерывности функции в
узлах, коэффициенты полинома выражают через вектор Ф
(е)
узловых
значений функции и координаты узлов и, проделав эквивалентные
преобразования, получают
eeе
ФN
(3)
где N
(е)
- матрица-строка, элементы которой называют функциями формы
конечного элемента.
Функции формы легко вычисляются в каждой точке конечного элемента
через координаты самой точки и координаты узлов элемента.
Этап 3. Объединение конечных элементов в ансамбль. На этом этапе уравнения
(4), относящиеся к отдельным элементам, объединяются в ансамбль, т. е. в
систему алгебраических уравнений:
NФ
(4)
Система (5) является моделью искомой непрерывной функции.
Этап 4. Определение вектора узловых значений функции. В общем случае вектор
Ф в (5) вначале неизвестен. Его определение - наиболее сложная процедура в
МКЭ.
Разработано несколько алгоритмов вычисления вектора Ф. Один из
алгоритмов основан на минимизации функционала, связанного с физическим
смыслом решаемой задачи, он состоит из 4-х этапов.
Найденные значения вектора Ф подставляют в (4), после чего значение
функции
легко вычисляется в любой точке заданной области.
Рассмотрим подробно этапы алгоритма МКЭ.
Этап 1: выделение конечных элементов.
Разбиение области на элементы - важный этап в МКЭ. От качества разбиения
во многом зависит точность получаемых результатов. Например, разбиение на
двухмерные элементы, близкие по форме к равносторонним треугольникам,
обеспечивает лучшие результаты по сравнению с разбиением на вытянутые по
форме треугольные элементы. Возможность легко варьировать размерами
элементов - важное свойство МКЭ.
Разбиение области на элементы обычно начинают от ее границы с целью
наиболее точной аппроксимации формы границы, затем производится разбиение
внутренних областей. Часто разбиение области на элементы производят в
несколько этапов. Сначала область разбивают на достаточно крупные подобласти
(подконструкции), границы между которыми проходят там, где изменяются
свойства материала, геометрия, приложенная нагрузка и пр. Затем каждая
подобласть взбивается на элементы. Резкого изменения размеров конечных
элементов на границах подобластей стараются избегать. На рис. 3. приведен
пример разбиения двухмерной области произвольной формы на треугольные
конечные элементы с криволинейными границами.
Рис. 3. Разбиение двухмерной области произвольной формы на треугольные конечные
элементы с криволинейными границами.
Нумерация узлов элементов (глобальная нумерация узлов) - следующая
процедура этапа выделения конечных элементов. Порядок нумерации имеет в
данном случае существенное значение, так как влияет на эффективность
последующих вычислений.
Матрица коэффициентов системы линейных алгебраических уравнений, к
которой приводит МКЭ, - сильно разреженная матрица ленточной структуры.
Ненулевые элементы такой матрицы располагаются параллельно главной
диагонали (рис. 4).
Рис. 4. Матрица ленточной структуры.
Целое число L, представляющее собой наибольшую разность между
номерами ненулевых элементов в строке, называется шириной полосы. Чем
меньше ширина полосы, тем меньший объем ОП требуется для хранения матрицы
при реализации МКЭ в САПР и тем меньше затраты машинного времени на
решение результирующей системы уравнений. Ширина полосы зависит, в свою
очередь, от числа степеней свободы узлов и способа нумерации последних.
Под числом степеней свободы понимают количество неизвестных функций,
определяемых в каждом узле. Так, например, для двухмерных задач гидравлики в
каждом узле определяют три переменные: давление и составляющие скорости по
осям X и Y.
При нумерации узлов предпочтителен способ, обеспечивающий
минимальную разность между номерами узлов в каждом отдельном элементе.
Если максимальную разность между номерами узлов для отдельного элемента
обозначить N, а число степеней свободы - М, то ширина полосы
MNL 1
В некоторых случаях уменьшение числа N может быть достигнуто
последовательной нумерацией узлов при движении в направлении наименьшего
размера рассматриваемой области.
На рис. 5 приведены два различных способа нумерации узлов произвольной
области, разбитой на конечные элементы. При первом способе (рис. 5, а) N=14,
при втором (рис. 5, б) N=5. Ширина полосы для представленных способов при
одной степени свободы в узле получается равной соответственно 15 и 6; при двух
степенях свободы - 30 и 12. Рациональная нумерация в случае рис. 5, б сокращает
необходимый объем ОП почти в три раза по сравнению со случаем рис. 5, б.
Рис. 5. Способы нумерации узлов при разбиении двухмерной области на конечные элементы.
Информация о способе разбиения области на конечные элементы и
нумерации узлов является исходной для всех следующих этапов алгоритмов МКЭ
при реализации метода в САПР. При этом требуется указывать не только номер,
но и координаты каждого узла и его принадлежность к определенным конечным
элементам. Такого рода информация называется топологической и обычно
содержит примерно в 6 раз больше чисел, чем количество узлов системы.
При описании области, разбитой на конечные элементы, необходимо
задавать: тип конечного элемента; его порядковый номер; номера узлов элемента;
координаты узлов, информацию о соединении элементов между собой; значение
физических параметров объекта в пределах каждого конечного элемента.
Промышленная эксплуатация программной системы долгое время тормозилась
именно сложностью подготовки исходных данных, объем которых в некоторых
случаях достигал нескольких сотен тысяч.
Усилия разработчиков программ МКЭ в составе САПР в последние годы
были направлены на создание подсистем автоматизированной подготовки
топологической информации, основу которых составляют специальные
программы, называемые препроцессорами. Препроцессоры либо непосредственно
включаются в состав программных комплексов, реализующих МК; либо
существуют в виде автономных программ.
Этап 2: определение аппроксимирующей функции элементов.
Эту процедуру можно выполнить один раз для типичного элемента области
безотносительно к его топологическому положению в ней. Полученная функция
используется далее для всех остальных элементов области того же вида. Эта
особенность является важным аспектом МКЭ. Благодаря ей элементы с однажды
определенными функциями легко включаются в библиотеку элементов
соответствующего программного комплекса. Далее эти элементы применяются
для решения разнообразных краевых задач.
В качестве аппроксимирующих функций элементов чаще всего используются
полиномы. В зависимости от степени полимома конечные элементы делятся на
симплекс-, комплекс- и мультиплекс-элементы. Полиномы симплекс-элементов
содержат константы и линейные члены; полиномы комплекс-элементов -
константы, линейные члены, а также члены более высоких степеней. Комплекс-
элементы, как правило, кроме граничных имеют дополнительные внутренние
узлы. Полиномы мультиплекс-элементов также содержат члены более высоких
степеней. На мультиплекс-элементы накладывается дополнительное условие: их
границы должны быть параллельны координатным осям.
Одномерный симплекс-элемент представляет собой отрезок (рис. 6).
Рис. 6. Одномерный симплекс-элемент.
При определении функции этого элемента для простоты будем считать, что
узловые значения искомой непрерывной функции, определенные на концах
отрезка, известны. По длине отрезка значение функции
аппроксимируется
полиномом:
xaa
21
(5)
Коэффициенты
1
a
и
2
a
определяются через узловые значения функции
i
и
j
в соответствии с условием непрерывности:
при
при
jj
ii
Xx
Xx
(6)
Подставив (6) в (5), получим систему уравнений:
jj
ii
XaaФ
XaaФ
21
21
;
решая, которую определим
1
a
и
2
a
:
L
ФФ
XXL
ФФ
XX
ФФ
aXaXaФФ
XXLXaФa
ij
ii
ij
ij
ij
jiij
ijii
222
21
;
; ;
L
XФXФ
L
XФФXXФ
L
XФФLФ
X
L
ФФ
Фa
ijjiiijijiiiji
i
ij
i
)()()(
1
т.е.
LФФaLXФXФa
ijijji
/ ;/
21
Подставив вычисленные значения коэффициентов аппроксимирующего
полинома в (5), получим
x
L
ФФ
L
XФXФ
ijijji
Проведем эквивалентные преобразования правой части:
L
Xx
Ф
L
xX
Ф
L
xФxФXФXФ
x
L
ФФ
L
XФXФ
i
j
j
i
ijijjiijijji
(7)
Члены полученного уравнения, заключенные в скобки, являются функциями
формы одномерного симплекс-элемента:
LxXN
ji
/
;
LXxN
ij
/
(8)
С учетом обозначений (8) уравнение (7) принимает вид
jjii
ФNФN
(9)
или в матричной форме
NФ
(10),
где
ji
NNN ,
- матрица-строка;
j
i
Ф
Ф
Ф
- вектор-столбец.
Функция формы обладает следующим свойством: функция формы с номером
i равна 1 в узле с соответствующим номером и равна 0 во всех других узлах. Не
представляет труда убедиться в наличии этого свойства у функций формы (8).
Этап 3: объединение конечных элементов в ансамбль.
Основу этого этапа составляет замена произвольно назначенных выше
номеров узлов i, j, k на номера, присвоенные узлам в процессе разбиения
рассматриваемой области. Эта процедура приводит к системе линейных
алгебраических уравнений, позволяющей при известных узловых значениях
искомой функции получить значение последней в любой точке области.
Рассмотрим процедуру составления ансамбля конечных элементов для
сформулированной выше задачи нахождения поля температур в стержне (рис. 2,
а). Кусочно-элементная модель области приведена на рис. 2, б, а функция
отдельного элемента определяется уравнением (10).
Можно написать следующее соответствие между произвольными номерами i,
j, фигурирующими в уравнении (10), и глобальными номерами узлов
рассматриваемой дискретной модели для
элемента 1:
2,1 ji
(11 а)
элемента 2:
3,2 ji
(11 б)
элемента 3:
4,3 ji
(11 в)
элемента 4:
5,4 ji
(11 г)
элемента 5:
6,5 ji
(11 д)
Подставив значения номеров узлов (11) в (9), получим:
6
5
65
5
5
5
5
4
54
4
4
4
4
3
43
3
3
3
3
2
32
2
2
2
2
1
21
1
1
1
TNTN
TNTN
TNTN
TNTN
TNTN
(12)
где верхние индексы в скобках относятся к номеру элемента.
В выражениях для функций формы элемента (9) значения произвольных
номеров i, j также следует изменить в соответствии с (11). Тогда значения
2
3
N
,
3
3
N
, например, определяются по формулам:
2
2
2
3
/ LXxN
3
4
3
3
/ LxXN
Очевидно, что
2
3
N
и
3
3
N
не равны друг другу даже в случае равенства длин
элементов L
(2)
и L
(3)
. При известных значениях узловых величин Т
1
-Т
6
уравнения
(12) позволяют определить значение температуры в любой точке стержня.
Расширенная форма описания моделей имеет некоторые преимущества при
реализации следующих этапов алгоритма МКЭ.
Определение вектора узловых значений функций. Для этой цели
используется несколько методов.
Метод, основанный на вариационной постановке задачи, требует
минимизации некоторого специально подобранного функционала, который связан
с физическим смыслом задачи. Подбор функционала является нетривиальной
процедурой, требующей глубоких знаний в конкретной предметной области.
Пример минимизации функционала в задаче о нахождении распределения температуры в
стержне (рис. 7).
Рис. 7. Пример минимизации функционала при нахождении распределения температуры в
стержне.
При указанном методе минимизируется функционал:
dSTTqTdV
dx
dT
F
SV
x
2
*
2
5,05,0
(13)
где V - объем тела; S - площадь границы.
В функционал F входят оба граничных условия (2). При минимизации функционала
используется множество функций элементов дискретизированной области. Для простоты
вычислений будем считать, что стержень разбит всего на два элемента (в практических случаях
этого недостаточно) (рис. 12). Тогда
2
1
21
1
1
1
TNTNT
;
3
2
32
2
2
2
TNTNT
(14)
Функционал (13) удобно представить в виде
2
2
*1
2
1 2
5,05,0 dSTTqTdSdV
dx
dT
F
S SV
x
(15)
где S
1
и S
2
- площади сечений стержня, на которых заданы граничные условия (2 а) и (2 б)
соответственно.
Для вычисления объемного интеграла в (15) его необходимо разбить на два слагаемых в
соответствии с принятой конечной элементной моделью:
21
2
2
2
1
2
1
2
5,05,05,0
V
x
V
x
V
x
dV
dx
dT
dV
dx
dT
dV
dx
dT
(16)
Производные в (16) вычисляются с учетом (14) и (8)
2
32
2
1
21
1
/
;/
LTT
dx
dT
LTT
dx
dT
(17)
Подставив (17) в (16) и считая, что
dxSdV
ee
получим
2
32
2
2
2
2
21
1
1
1
2
5,05,05,0 TT
L
S
TT
L
S
dV
dx
dT
xx
V
x
Второе и третье слагаемые в (15) вычисляются просто так как подынтегральным
функциям соответствуют узловые значения Т
1
и Т
3
:
1
111
S
SqTqTdS
;
2
*3*
2
322
2
*
25,05,0
2
TTTTSdSTT
S
где S
1
и S
2
- площади поверхностей, на которых заданы q и
(для рассматриваемого примера
1
1
SS
и
2
2
SS
).
Значение функционала F вычисляется простым суммированием последних трех
выражений:
2
*3*
2
3211
2
332
2
2
22
221
2
1
1
25,0
25,025,0
TTTTSTqS
TTTTCTTTTCF
(18)
где
1
1
11
/ LSC
x
и
2
2
22
/ LSC
x
Для минимизации функционала F необходимо выполнение условий:
0][/
;0][/
;0/
*232
)2(
2
)2(
3
3
)2(
2
)2()1(
1
)1(
2
12
)1(
1
)1(
1
TSTSCTCTF
TCTCCTCTF
qSTCTCTF
или в матричной форме
*2
1
3
2
1
2
)2()2(
)2()2()1()1(
)1()1(
0
0
0
TS
qS
T
T
T
SCC
CCCC
CC
(19)
В общем виде (19) можно представить так:
ВКТ
, что соответствует (7).
Примечание. Матрица коэффициентов К в (19) называется матрицей жесткости, хотя
по физическому смыслу данной задачи ее удобнее было бы назвать матрицей
теплопроводности. Такое название матрицы К пришло из строительной механики, где МКЭ
начал применяться раньше, чем в других областях техники.
Зная характеристики материала, из системы (19) можно определить узловые значения T
1
,
Т
2
, Т
3
. Из (18) и (19) нетрудно заметить, что однотипные конечные элементы вносят в эти
выражения слагаемые одного вида. Поэтому при реализации МКЭ в САПР вклад элемента
определенного типа в матрицу жесткости вычисляется только один раз, а затем используется во
всех необходимых случаях.
Метод Галеркина - другой широко известный метод вычисления вектора
узловых значений - представляет собой частный случай более общего метода
взвешенных невязок. Основным преимуществом этого метода является то, что
основой для него служит исходное дифференциальное уравнение. Поэтому метод
Галеркина с успехом применяется при решении задач, из которых не удается
подобрать функционал для минимизации (например задач, математическим
описанием которых служат уравнения Навье-Стокса).
Метод Галеркина основан на минимизации ошибки
fLu
приближенного
решения u исходного дифференциального уравнения
0 fL
, где L -
дифференциальный оператор.
Для минимизации
в заданной области G требуется выполнение равенства
G
i
dGN 0
для каждой из функций N
i
.
Сочетание метода Галеркина с МКЭ приводит к системе уравнений
G
kjidGLN ..., , , , при 0
где
L
- левая часть исходного дифференциального уравнения, описывающего
непрерывную функцию
. Технику получения разрешающей системы уравнений
методом Галеркина легко проиллюстрировать на примере задачи об отыскании
температурного поля в однородном стержне (рис. 1).
Применив метод Галеркина к (1), получим
V
x
t
dV
dx
Td
N 0
2
2
(20)
Подставим в (33) формулу дифференцирования произведений:
V
x
t
V V
x
t
x
t
dV
dx
dT
dx
dN
dV
dx
dT
N
dx
d
dV
dx
Td
N
2
2
(21)
Интерполяционная функция Т не сохраняет постоянства по длине стержня,
поэтому интеграл в (21) можно представить суммой соответствующих интегралов
для отдельных элементов. Так, второй интеграл в (21) можно представить в виде
2
1e
V
e
e
e
te
V
x
t
dV
dx
dT
dx
dN
dV
dx
dT
dx
dN
x
Вычислим в этом уравнении интегралы, относящиеся к отдельным
элементам:
1
1
1
e
e
i
j
te
L
L
Xx
xX
dx
d
dx
dN
(22)
j
i
e
ee
e
T
T
L
TN
dx
d
dx
dT
11
1
(23)
С учетом (22) и (23)
j
i
e
ee
V
e
j
i
ee
e
V
e
e
e
t
e
T
T
L
S
dV
T
T
LL
dV
dx
dT
dx
dN
xx
x
11
11
11
11
(24)
Первый интеграл в (21) на основании теоремы Остроградского-Гаусса
преобразуется к виду
dSl
dx
dT
NdV
dx
dT
N
dx
d
x
V S
x
t
x
t
(25)
где
dndTdxdTl
x
//
; n - внешняя нормаль к рассматриваемой поверхности.
С учетом краевого условия (2 а) в точке для первого элемента интеграл (25)
принимает вид
dSq
L
X
L
X
dS
dn
dT
N
S S
x
1
1
1
2
1
1
(26)
Подставив значения
0
i
X
,
1
2
LX
в (26), получим
00
1
1
1
qS
dS
q
dS
dn
dT
N
S S
x
t
(27)
С учетом краевого условия (2 б) в точке
2
Lx
второго элемента интеграл
(25) запишется так:
*32*3
2
22
1
0
1
TTSdSTT
q
dS
dx
dT
N
SS
x
t
(28)
Просуммировав выражения вида (24) для первого и второго элементов и
выражения (27) и (28) и приравняв сумму нулю, получим систему уравнений
*2
1
3
2
1
2
22
2211
11
0
0
0
TS
qS
T
T
T
SCC
CCCC
CC
(29)
где
22221111
/ и / LSCLSC
xx
Завершающим шагом этапа определения вектора силовых значений Ф
является решение системы линейных алгебраических уравнений.
Примечание. Основные особенности этого шага - большая размерность и сильная
разреженность матрицы коэффициентов системы. В связи с этим для реализации МКЭ в САПР
разработаны специальные способы хранения матрицы жесткости, позволяющие уменьшить
необходимый для этого объем ОП. Для нахождения узловых значений функций применяются
методы преобразования и решения системы, направленные на снижение затрат машинного
времени.
Метод граничных элементов.
При решении краевых задач приближенные модели технических объектов
можно строить на основе интегральных уравнений. При этом первый шаг на пути
к решению состоит в переходе от дифференциальных уравнений в частных
производных к эквивалентным интегральным уравнениям. Во многих случаях,
когда такой переход оказывается успешным, решение исходной задачи может
быть получено с минимальными вычислительными затратами и высокой
степенью точности. Кроме того, размерность исходной задачи понижается на 1,
двухмерные задачи преобразуются в одномерные.
Примером указанного подхода к решению краевых задач служат методы
интегральных граничных элементов (МГЭ). Развитие МГЭ началось сравнительно
недавно, причем мощным стимулом к этому послужило создание
быстродействующих ЭВМ. Все разновидности МГЭ используют принцип
суперпозиции, поэтому область их применения ограничена классом полностью
линейных или линейных относительно приращений задач. Однако к такому
классу относятся многие важные для развития техники задачи, например МГЭ
успешно используются для решения задач теории упругости, механики жидкости
и газов.
Все МГЭ строятся на основе общих принципов. При этом различают прямые
и непрямые МГЭ.
В прямых МГЭ искомыми переменными краевой задачи являются величины,
имеющие реальный физический смысл, например, в задачах теории упругости -
усилия и смещения, возникающие в элементах конструкции.
В непрямых МГЭ решение исходной задачи выражается через функции
плотности, которые сами по себе не имеют реального физического смысла. После
того как функции плотности найдены, значения реальных физических параметров
задачи могут быть получены из них путем простого интегрирования.
В любом варианте МГЭ результатом перехода от дифференциальных
уравнений в частных производных к интегральным уравнениям в конечном счете
является система уравнений, включающая значения переменных только на
границе заданной области. Поэтому в отличие МКЭ и МКР последующая
дискретизация задачи проводится только на границе исследуемой области.
Последнее обусловливает, во-первых, более высокую по сравнению с МКР и
МКЭ точность решения, во-вторых, существенно меньший объем входных
данных при реализации методов на ЭВМ.
Переход от исходного дифференциального уравнения интегральному.
Рассмотрим на простом примере алгоритм перехода. В двухмерной
однородной области произвольной формы с коэффициентом проницаемости k
требуется найти распределение функции
, описанной уравнением
0//
2
2
22
1
2
xx
(1),
которое является частным случаем квазигармонического уравнения
0
Q
z
K
zy
K
yx
K
x
zyx
(1.1) при
,
21
yххх
.
На границе L рассматриваемой области заданны граничные условия первого
рода
00
xgx
(2)
Этап 1. Нахождение сингулярного решения. В МГЭ используется то
обстоятельство, что для большинства уравнений в частных производных