лекция. глава4_5_6_2018. Макро и микромоделирование глава 4 Метод конечных элементов

• Название: Макро и микромоделирование глава 4 Метод конечных элементов
• Анкор: лекция
• Дата: 16.04.2023
• Размер: 1.17 Mb.
• Имя файла: глава4_5_6_2018.pdf
• Тип: Документы; #1065236
• страница: 1 из 5

114
РАЗДЕЛ 2 МАКРО- И МИКРОМОДЕЛИРОВАНИЕ
Глава 4 Метод конечных элементов
В настоящее время МКЭ является одним из наиболее популярных ме- тодов решения краевых задач в САПР. В математическом отношении МКЭ относится к группе вариационно-разностных. Строгое доказательство та- ких важных свойств, как устойчивость, сходимость и точность метода, проводится в соответствующих разделах математики и часто представляет собой непростую проблему. Тем не менее, метод получил широкое распро- странение для решения задач механики сплошных сред, благодаря своей универсальности, ясной инженерной формализации и удобству реализации на ЭВМ [5-9].
МКЭ был впервые применен в инженерных приложениях в начале
1950-х годов. В 1956г. группа Тернера из Boeing Aircraft Co описала про- цедуру, включающую некоторые характерные черты МКЭ. Интересно, что
МКЭ был развит независимо в прикладной математике. В 1943г. Курант описал процедуру решения, основанную на принципе минимума потенци- альной энергии. Но до 1960-х годов математики и инженеры не предпри- нимали совместных попыток совершенствования этого метода.
Распространение МКЭ на другие задачи было предпринято в начале
1960-х годов на основе вариационного подхода. Буквально в последние три десятилетия начали применять другие формулировки МКЭ. Наиболее популярные из них - метод Галеркина, метод наименьших квадратов, пря- мой метод, метод глобального баланса или метод Одена.
В принципе, несмотря на большое разнообразие в формулировках,
МКЭ может быть охарактеризован следующими этапами:
1. Физическая область задачи делится на подобласти или конечные элементы (этап дискретизации).
2. Зависимая переменная аппроксимируется функцией специального вида на каждом конечном элементе и, следовательно, по всей области. В этом случае параметры этих аппроксимаций в дальнейшем считаются ос- новными неизвестными задачи (этап аппроксимации).

115 3. Подстановка аппроксимирующих коэффициентов в определяющие уравнения дает систему линейных алгебраических уравнений (ЛАУ), ре- шая которые можно найти основные неизвестные и, следовательно, полу- чить приближенное решение задачи (этап алгебраизации).
4.1 Общая вариационная формулировка МКЭ
Пусть поведение искомой функции

(x,y,z) внутри заданной области V с границей S описывается некоторым дифференциальным уравнением 2m - го порядка:
)
,
,
(
)
,
,
,
,
(
2
z
y
x
Q
z
y
x
k
L
m


,
(4.1) где К – параметр, обобщенно характеризующий свойства сплошной среды;
Q – внешнее воздействие; L
2m
– дифференциальный оператор порядка 2m.
Уравнение (4.1) дополняется совокупностью m -краевых условий на поверхности S, ограничивающей область V.
При решении задач механики сплошных сред для определения функ- ции

(x,y,z) вместо рассмотрения уравнения (4.1) часто используют вариа- ционный подход. Для рассматриваемого класса задач можно показать, что решение

(x,y,z) уравнения (4.1) совпадает с функцией, минимизирующей некоторый функционал F(

) , который содержит производные от

(x,y,z) до m -го порядка. Это облегчает подбор аппроксимирующих функций для

(x,y,z), поскольку для получения значения функционала F(

) требуется обеспечение непрерывности функции

(x,y,z) и ее производных лишь до
(m-1) -го порядка включительно.
Такая трактовка обусловливает следующую последовательность прове- дения расчета по МКЭ:
I. Дискретизация области решения. Разбиение области решения на ко- нечные элементы является первым шагом на пути к решению задачи. Этот шаг не имеет теоретического обоснования. Искусство разбиения области зависит от имеющихся инженерных навыков и опыта моделирования. Сле- дует отметить, что несовершенная дискретная модель будет приводить к

116 значительным погрешностям расчета, если даже все остальные этапы ме- тода выполнены с достаточной точностью.
Выбор типа, формы элемента и числа его узловых точек зависит от ха- рактера рассматриваемой задачи и от той точности решения, которую тре- буется обеспечить. Рекомендации к обоснованному использованию конеч- ных элементов будут сформулированы в дальнейшем. Здесь же следует отметить, что при замене исходной конструкции его дискретной моделью стараются обеспечить как можно большую идентичность в поведении кон- струкции и ее модели.
2. Аппроксимация искомой функции. В качестве основных неизвест- ных в МКЭ принимаются узловые значения искомой функции и ее частные производные до m -го порядка. Например, при решении краевых задач, описываемых квазигармоническим уравнением (2m=2), в качестве неиз- вестных в каждой i-той узловой точке достаточно принять значение опре- деляемой функции φ
i
Общее число неизвестных определяет число степеней свободы, от ко- торого зависит точность определения искомой функции в объеме каждого конечного элемента, а следовательно, и во всей области решения задачи.
После выбора узловых неизвестных строится аппроксимирующий по- лином, который выражает закон изменения искомой функции φ(x,y,z) по объему конечного элемента через значения его узловых неизвестных. По- лученные полиномы должны обеспечить непрерывность функции φ(x,y,z) и eе производных до (m-1) -го порядка включительно во всей области реше- ния. В каждом из полиномов должны содержаться члены, обеспечивающие их переход к постоянным значениям при уменьшении размеров конечного элемента.
В дальнейшем вопросы построения аппроксимирующих полиномов для типичных конечных элементов будут рассмотрены подробно. Пока же примем, что значение непрерывной функции φ
(е)
в произвольной точке е-го конечного элемента аппроксимируется полиномом.
)
(
)
(
)
(
}
{
]
[
e
e
e
N



,
(4.2)

117 где [N]
(e)
- матрица-строка, элементы которой называются функциями формы конечного элемента; {

}
(e)
- вектор узловых неизвестных e -го ко- нечного элемента. Тогда аппроксимация закона изменения искомой функ- ции φ(x,y,z) по всей области решения определяется суммой




M
e
e
e
N
N
z
y
x
1
)
(
)
(
}
]{
[
}
{
]
[
)
,
,
(



,
(4.3) где
М
- число конечных элементов дискретной модели;
[N]=[[N]
(1)
[N]
(2)
..[N]
(N)
]; {

}- вектор узловых неизвестных для всей сово- купности конечных элементов.
Система (4.3) является моделью искомой непрерывной функции.
3. Получение основной системы разрешающих уравнений.
В общем случае вектор {

} неизвестен. В вариационном МКЭ алго- ритм получения основной системы разрешающих уравнений основан на минимизации функционала F(φ) и состоит из следующих этапов:
Этап I. Выбор функционала F , который представляется суммой соот- ветствующих функционалов, относящихся к отдельным конечным элемен- там



M
e
e
F
F
1
)
(
,
(4.4) где F
(e)
- элементный вклад в функционал F.
Этап 2. Подстановка аппроксимирующего выражения (4.2) в уравнение функционала F(φ) и получение элементных вкладов в функционал F.
Этап 3. Минимизация по вектору {

} функционала F. Для этого со- ставляются уравнения
0
}
{
)
(



e
e
F

(4.5)
Если справедливо выражение (4.4), то суммирование выражений (4.5) по конечным элементам приводит к системе алгебраических уравнений

118
}
{
}
]{
[
R
K


,
(4.6) где [К] - матрица жесткости системы элементов; {R}- вектор нагрузки.
Этап 4. Решение системы (4.6), позволяющее определить неизвестный вектор узловых значений. Для линейных краевых задач система уравнений
(4.13) линейна. Для ее решения oбычно используются методы Гаусса, Хо- лесского, сопряженных градиентов и иногда, при очень высоком порядке системы, итерационные методы.
Для нелинейных краевых задач система (4.6) нелинейна, поскольку матрица [К] является функцией определяемых неизвестных параметров
{

}. При решении нелинейной системы алгебраических уравнений ис- пользуются итерационные методы.
Определение выходных параметров краевой задачи. Пусть вектор {

} найден. Тогда с помощью зависимости (4.2) можно определить вектор
{

}
е
, для каждого конечного элемента. Компоненты напряженно- деформированного состояния, которые настакже могут интересовать при решении краевых задач, определяются либо дифференцированием полу- ченного выражения для φ
(е)
либо непосредственно через узловые значения искомых npoизводных, если последние входят в состав вектора {

}.
4.2 Вариационная формулировка двумерного МКЭ
Рассмотрим двумерную задачу теплопроводности через брус квадрат- ного сечения (рисунок 4.1).
Требуется найти распределение температуры в брусе и в частности температуру в точке А.

119
T=100
C
L=2
A
L=1
L=1
/4
L=2
T=50 C

T/

x
.
0
Рисунок 4.1 - Брус квадратного сечения
Воспользуемся вариационной формулировкой МКЭ, в которой приме- няется глобальная система координат
ОХУ
В данном случае определяющим уравнением является уравнение
Лапласа
0 2
2 2
2






y
T
x
T
,
(4.7) с граничными условиями Дирихле на части границы
,
100
,
50




T
T
,
0
L
y
y


(4.8) и условиями Неймана на остальной части границы
,
0
,
0






x
T
x
T
,
0
L
x
x


(4.9)
Задача определена на множестве R, состоящем из D и границы S, т.е.
R=D+S.
Уравнения (4.7) - (4.9) не будут использоваться, непосредственно вме- сто них будет построена эквивалентная вариационная формулировка. С помощью вариационного исчисления можно показать [7], что решение

120
Т(х,у), удовлетворяющее уравнениям (4.7) - (4.9), совпадает с функцией, которая минимизирует функционал
 








D
dxdy
y
T
x
T
]
)
(
)
[(
2 1
2 2

,
(4.10) где

T
(x, y) - функция из доступного множества пробных функций, задан- ных в области D.
Для этой задачи пробные функции

T (x,y) являются допустимыми, если они непрерывны и имеют кусочно-непрерывные первые производные.
Кроме того, пробные функции должны удовлетворять главным граничным условиям (4.8). Граничные условия Неймана (4.9) будут выполняться ав- томатически для функции, минимизирующей функционал (4.10) как есте- ственное следствие вариационной формулировки и будут называться есте- ственными граничными условиями.
Следует заметить, что для большинства задач, возникающих в процессе проектирования, вариационная формулировка известна.
Решение задачи состоит из трех этапов.
Этап I. Разобьем область на l конечных элементов. Обычно это разбие- ние называется триангуляцией области, так как простые элементы часто являются треугольниками. Общее число узлов обозначим n (рисунок 4.2).
L=2
L=2
0
l
1
l
2
l
4
l
6
l
3
l
5
l
1
l
i
l
i-1
l
i-2
Рисунок 4.2 - Разбиение области на конечные элементы

121
Разбиение области и условия непрерывности, накладываемые на проб- ные функции, позволяют записать функционал (4.10) в виде



l
i
ei
1


,
(4.11) где

ei
– элементарный вклад, определяемый равенством
 


























ei
ei
ei
ei
dxdy
y
T
x
T
2 2
2 1

,
(4.12)
Для определения элементного вклада рассмотрим типичный элемент e
i
(рисунок 4.3), который характеризуется координатами узлов с номерами узлов.
Этап 2. Построение пространства пробных функций.
Для произвольного элемента e
i
в этом примере пробная функция T
ei
(x,y) выбирается линейной, т.е.
y
x
y
x
T
ei
3 2
1
)
,
(






, i
e y
,
x

(4.13) где

1
,

2
,

3
- постоянные, в общем случае различные для разных элемен- тов.
С целью определения этих постоянных запишем последовательно уравнение (4.13) для узлов i, j, m
,
,
3 2
1 3
2 1
3 2
1
m
m
m
j
j
j
i
i
i
y
x
T
y
x
T
y
x
T


















(4.14) где T
i
, T
j
, T
m
– значение T в узлах i, j, m.
Система уравнений (4.14) имеет единственное решение, т.к. определи- тель ее коэффициентов не равен нулю, т.е.

122 0
1 1
1 2













m
m
j
j
i
i
y
x
y
x
y
x
(4.15)
Этот определитель равен удвоенной площади треугольника, т.к. пло- щадь треугольника никогда не равна нулю, то решение

1
,

2
и

3
суще- ствует и единично. Решая (4.14), получим для

1
,

2
и

3
следующие выра- жения:
),
)(
2 1
(
),
)(
2 1
(
),
)(
2 1
(
3 2
1
m
m
j
j
i
i
m
m
j
j
i
i
m
m
j
j
i
i
T
c
T
c
T
c
T
b
T
b
T
b
T
a
T
a
T
a















(4.16) где
j
m
m
ji
i
y
x
y
x
a


;
m
j
i
y
y
b


;
j
m
i
x
x
c


,
(4.17) а постоянные a
j
, a
m
, b
j
, b
m
, c
j
, c
m
могут быть получены циклической пере- становкой индексов. Подстановка выражений (4.16) -в (4.13) дает следую- щие представления через базисные функции:
]
)
(
)
(
)
[(
2 1
)
,
(
m
m
m
m
j
j
j
j
i
i
i
i
ei
T
y
c
x
b
a
T
y
c
x
b
a
T
y
c
x
b
a
y
x
T










(4.18) или
e
e
m
m
j
j
i
i
ei
T
N
T
N
T
N
T
N
y
x
T
}
{
]
[
)
,
(




,
(4.19) где
]
[
]
[
m
j
i
e
N
N
N
N

- является матрицей базисных функций или функций формы;
}
T
T
T
{
}
T
{
m j
i

- вектор узловых значений. Получим производные, входящие в функционал (4.12),
).
(
2 1
),
(
2 1
m
m
j
j
i
i
ei
m
m
j
j
i
i
ei
T
c
T
c
T
c
y
T
T
b
T
b
T
b
x
T












(4.20)
Подстановка (4.19) в выражении для элементного вклада (4.12) дает

123
 







ei
m
m
j
j
i
i
m
m
j
j
i
i
ei
dxdy
T
c
T
c
T
c
T
b
T
b
T
b
x
]
)
(
)
[(
8 1
2 2
2
. (4.21)
Выражение под интегралом не зависит от x и y



ei
dxdy
,
(4.22) тогда равенство (4.21) переписывается следующим образом:
]
)
(
)
[(
8 1
2 2
2
m
m
j
j
i
i
m
m
j
j
i
i
ei
T
c
T
c
T
c
T
b
T
b
T
b








(4.23)
Выражение (4.23) может быть получено для элемента. Подставляя все эти элементные вклады в (4.11), преобразуем функционал, заданный ра- венством (4.10), в функционал всех узловых значений
T
1
, T
2
,

, T
n
,т.е.

=

(T
1
, T
2
,

, T
n
).
(4.24)
Здесь узловые параметры рассматриваются в качестве переменных, значение которых необходимо определить. Условия минимума могут быть записаны в виде
0











p
T

, p=1,2,…,n.
(4.25)
Подстановка (4.11) в (4.25) позволяет представить эти уравнения сле- дующим образом:








l
i
p
ei
p
T
x
T
x
1 0
, p=1,2,…,n.
(4.26)
Очевидно, что при суммировании в (4.26) ненулевой вклад дают только те элементы, которые содержат узел р. Если узлы элемента (i, j и m) имеют номера в глобальной системе p, q и r соответственно, то дифференцирова-

124 ние (4.12) по Т
p
приводит к выражению
ei
ei
p
r
r
q
q
p
p
p
r
r
q
q
p
p
p
p
ei
T
K
T
c
T
c
T
c
c
T
b
T
b
T
b
b
T
x
}
{
}
{
)]
(
)
(
[
4 1










.(4.27)
Объединение компонент элементных уравнений, задаваемое равен- ством (4.26), называется объединением по узлам. Это объединение позво- ляет получить систему уравнений, решение которой дают нам узловые значения температуры.
Выражение (4.26) определяют сущность 3 этапа МКЭ.
Рассмотрим процедуру поэлементного объединения. Определим вклад е -го элемента в общее системное уравнение
e
e
e
m
j
i
e
mm
mj
mi
jm
jj
ji
im
ij
ii
m
e
j
e
i
e
e
e
T
K
T
T
T
K
K
K
K
K
K
K
K
K
T
T
T
T
}
{
]
[
}
{





















































,
(4.28) где [K]
e
- матрица жесткости элемента.
Процедуру объединения в этом случае можно записать как






l
e
e
e
T
K
T
1 0
}
{
]
[
}
{

или
0
}
]{
[
}
{




T
K
T
x
,
(4.29) где [К] - матрица жесткости системы элементов, {Т}
T
= {T
1
, T
2
,…, T
n
}.
Уравнение (4.29) может быть решено последовательным исключением либо одним из стандартных методов решения систем ЛАУ, после учета граничных условий (4.8).
Таким образом, вариационная формулировка МКЭ позволяет получить системное уравнение МКЭ (4.29). Матрица жесткости [К] - типичная раз- реженная матрица ленточной структуры.
Подготовка к решению задачи МКЭ требует следующее:

125 1) разбить область и на треугольники (в общем случае это может быть любой многоугольник);
2) перенумеровать вершины таким образом, чтобы матрица жесткости имела ленточную структуру. Ширина ленты определяется наибольшим разрывом между номерами соседних вершин.
Для случаев, когда узлу соответствует более чем один параметр, как, например, у эрмитовых элементов полуширина матрицы [K] определяется как
1
)
1
(



q
b
S

, гдеα - максимальная разность между самым мень- шим и самым большим номерами узла любого элемента системы;
3) вычислить соответствующие интегралы (4.21), что довольно просто, если использовать нормализованные (местные) координаты в пределах каждого треугольника элемента.
Заметим следующее: а) простой элемент является треугольником К;б) число степеней свободы определяется тремя вершинами треугольника;, в) пространством интерполирующих функций (пробных) является простран- ство полиномов Р
1
степени ≤ 1. Триплет {K, Σ K, P
1
} /5/ называется тре- угольным конечным элементом.
4.3 Метод конечных элементов в задачах теории упругости
Рассмотрим тело, определенное двумерной областью S и границей C с заданными граничными условиями (рисунок 4.3). При действии внешних сил для данного случая можно воспользоваться основным уравнением тео- рии упругости.
1. Уравнение равновесия (Лавье)





















,
0
,
0
Y
y
x
X
y
x
y
xy
xy
x




(4.30) где X и Y – объемные силы;
 


T
xy
y
x





- вектор столбец перемеще- ний.
2. Зависимости между деформациями и перемещениями

126
(уравнение Коши)
x
U
x




,
y
V
y




,
y
V
x
U
xy







,
(4.31) где U и V– соответственно перемещения точки в направлении x и y;
 


T
xy
y
x





- вектор-столбец деформаций. а)
б) а) тело, определенное двумерной областью S и границей C с заданны- ми граничными условиями; б) напряженное состояние элемента dxdy
Рисунок 4.3 - Пример плоской задачи теории упругости
3. Уравнение состояния (закон Гука)














,
),
(
1
),
(
1
;
'
;
'











xy
xy
x
y
y
y
x
x
E
E
(4.32) где '
E
и
;

- соответственно модуль упругости при растяжении и коэф- фициент Пуассона.
Для плоского напряженного состояния
E
E

'
и



'
, для плоской деформации
)
1
/(
2
'



E
E
и
)
1
/(
'






127 4.
Граничные условия:
- геометрические, на С
u
:
U
U

,
V
V

,
- силовые, на С
σ
:


X
X

,


Y
Y

, где С
u и С
σ
– границы, на которых заданы условия;

X
и

Y
- поверхностные усилия;
V
- вектор, направленный по нормали к границе С
σ
наружу.
Эти силы можно представить в виде:





sin cos
xy
x
X


,





sin cos
y
xy
Y


(4.33)
Для рассматриваемой задачи можно показать, что ее решение совпада- ет с функцией, минимизирующей функционал:
















s
C
s
xy
xy
y
y
x
x
tdC
V
Y
U
X
tdxdy
V
Y
U
X
tdxdy
F













2 1
,
(4.34) где t – толщина пластины;
dC – элемент длины на границе C.
В качестве основных неизвестных приняты функции перемещения
U(x,y) и V(x,y), поэтому данная формализация носит название метода пере- мещений.
Рассмотрим основные этапы МКЭ для решения плоской задачи теории упругости на примере использования треугольного симплекс-элемента
(рисунок 4.4).
Рисунок 4.4 - Аппроксимация области конечными элементами

128
Для аппроксимации искомой функции воспользуемся выражением:
   
 
 
e
e
T
e
N
UV
f



,
 


T
k
k
j
j
i
i
e
V
U
V
U
V
U


, (4.35) где
 
e

- вектор узловых перемещений е-го элемента;
 
e
N
- матрица функций формы.
Выбор элементов и функций формы оказывает существенное влияние на точность решения.
Если выражение (4.35) подставить в зависимость (4.31), связывающие перемещения и деформации, и провести соответствующие преобразования, то
 
 
 
e
e
B



,
(4.36) где
 
B
- матрица градиентов.
Для упругого тела воспользуемся законом Гука (4.32), который в мат- ричном виде можно представить:
 
 
   
)
(
0
e
e
e
D





,
(4.37) где
 
D - матрица упругости;
 
e
0

- вектор-столбец начальных деформа- ций.
Можно показать, что
 















2
)
1
(
0 0
0 1
0 1
1
'
'
'
2
'
'




E
D
(4.38)
Если зависимости (4.36) и (4.37) подставить в функционал (4.34) и применить процедуру его минимизации, то получим:

129
 
 
 
 
 
e
p
e
e
e
e
R
R
K
R



0


,
(4.39) где
 
    


e
S
T
e
tdxdy
B
D
B
K
,
(4.40)
 
   
 



e
S
T
e
tdxdy
D
B
R
0 0


,
(4.41)
 
 
 



e
S
T
e
p
tdxdy
p
N
R
,
 
 
T
Y
X
p

,
(4.42) где
 
e
R
- вектор эквивалентных узловых сил.
Уравнения равновесия для системы конечных элементов определяются в виде:
 
   
 
 
p
R
R
R
K



0



,
(4.43) где
 

R
- вектор внешних нагрузок.
Решая эту систему линейных алгебраических уравнений при заданных граничных условиях, можно определить перемещения в узлах и, следова- тельно, по формулам (4.35) - (4.37) перемещения, деформации и напряже- ния внутри элементов.
Следует отметить, что при решении трехмерных задач структура ос- новных матричных соотношений МКЭ не поменяется и их решение не вы- зывает принципиальных затруднений.
В данной работе в качестве базовых элементов используются треуголь- ные элементы при решении плоских задач и тетраэдральные элементы для решения трехмерных задач. При решении ряда задач применяются изопа- раметрические элементы первого и второго порядков (рисунки 4.5, 4.6).

130