Файл: Макро и микромоделирование глава 4 Метод конечных элементов.pdf

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

131
1
8
X
Y
0
5
5
2
3
4
1
Рисунок 4.6 - Изопараметрический элемент второго порядка
Треугольный симплекс-элемент /5, 7/. Допустим, что перемещения то- чек внутри элемента выражаются полиномом первого порядка от их коор- динат x и y.









y
x
v
y
x
u
6 5
4 3
2 1






(4.45)
Тогда выражение для перемещений произвольной точки внутри эле- мента имеет вид
e
k
j
i
N
I
N
I
N
I
f
}
]{
[
}
{





,
(4.46) где I – единичная матрица размерности 2х2, а
S
y
c
x
b
a
N
i
i
i
i
2




(4.47) и так далее.
Здесь S – площадь элемента,
j
k
k
i
i
y
x
y
x
a


,
k
i
i
y
y
b


,
j
k
i
x
x
c



132
Выбранная функция перемещений автоматически гарантирует непре- рывность перемещений между смежными элементами.
Изопараметрический квадратичный элемент [5,6]. Введем в изопара- метрическом четырехугольном элементе (рисунок 4.7) локальную коорди- натную систему ξ, η, которая удовлетворяет условиям
1 1




,
1 1




Функциям формы для интерполяции перемещений по их узловым зна- чениям на основе [5,6] имеют вид:
).
1
)(
1
(
2 1
),
1
)(
1
(
2 1
),
1
)(
1
(
4 1
)
1
)(
1
(
4 1
)
1
)(
1
(
2 1
2 0
0 2
2 0
0 2
0 0

























i
i
i
N
N
N
(4.48)
Здесь i – номер узла элемента (i = 1,2,…,8), а
i



0
,
i



0
,
ξ
i
, η
i
– значения локальных координат в i –м узле.
(X
1
,Y
1
)
X
Y
0
4
e
1
(X
2,
Y
2
)
2
3
(X
3,
Y
3
)
(X
4,
Y
4
)
Рисунок 4.7 - Четырехугольный элемент
Тетраэдральный элемент [6] . Перемещения любой точки определяются тремя компонентами u, v, w (рисунок 4.7). Тогда вектор перемещений име- ет вид
T
z
y
x
w
z
y
x
v
z
y
x
u
f
)}
,
,
(
)
,
,
(
)
,
,
(
{
}
{


133
По аналогии с (4.31) можно записать, например,
z
y
x
u
4 3
2 1








(4.49)
Перемещения произвольной точки в объеме элемента можно предста- вить в виде
e
m
k
j
i
N
I
N
I
N
I
N
I
f
}
]{
[
}
{






,
(4.50) где I – единичная матрица размерности 3х3, а
v
z
d
y
c
x
b
a
N
i
i
i
i
i
6





(4.51) и так далее.
Здесь V – объем конечного элемента; а коэффициенты a
i
, b
i
, c
i
и d
i
определяются через значения координат узловых точек элемента.
Трехмерный изопараметрический элемент с линейным полем переме- щений [5,6]. Удобнее всего выразить функции формы через координаты узлов на границе элемента и использовать локальные координаты. Они по- казаны на рисункe 1.6 и выбраны так, чтобы стороны конечного элемента совпадали с координатными линиями ±1.
Тогда получим следующие функции формы:
,
}
]{
8 2
1
[
}
{
e
N
I
N
I
N
I
f





(4.52) где I – единичная матрица размерности 3х3, а
)
1
)(
1
)(
1
(
8 1
0 0
0







i
N
,
i



0
,
i



0
,
i



0
,
(4.53)

134 где ξ
i
, η
i
, ζ
i
– значения локальных координат в i – м узле (i = 1, … ,8).
Отметим, что объемный элемент с 24-мя cтепенями свободы позволяет существенно увеличить число неизвестных параметров в выражении для компонентов перемещений (4.51) и тем самым обеспечить лучшую ап- проксимацию действительного НДС в пределах объема конечного элемен- та.
Описанные интерполирующие полиномы для каждого конечного эле- мента обеспечивают непрерывность функции U(x,y,z) и ее производных до
(m-1) - го порядка включительно. Здесь 2m - порядок определяющего диф- ференциального уравнения.
Опишем процедуры построения матриц жесткости рассматриваемых.
Треугольный симплекс-элемент. Для того чтобы записать закон Гука для плоской задачи в матричной форме, представим матрицу упругости в виде













2
)
1
(
0 0
0 1
0 1
1
]
[
2




E
D
(4.54)
Такая запись подразумевает, что деформации объединены в вектор в следующем порядке:
}
{
}
{
xy
y
x





Матрица [B] получается при помощи дифференцирования соотноше- ний (4.45) и ее можно представить в виде
]
[
]
[
k
j
i
B
B
B
B

(4.55) или











k
k
j
j
i
i
k
j
i
k
j
i
B
C
B
C
B
C
C
C
C
B
B
B
S
B
0 0
0 0
0 0
2 1
]
[
(4.56)