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

135
Отметим, что в этом случае матрица [B] не зависит от координат точек внутри элемента и, следовательно, деформациив нем постоянны.
Подставив теперь матрицы [B] и [D] в основное уравнение матрицы жесткости и узловых сил и проинтегрировав по объему конечного элемен- та, получим формулу для подсчета матрицы жесткости треугольного эле- мента
st
B
D
B
K
T
e
]
][
[
]
[
]
[

,
(4.57) где S - площадь элемента; t - толщина элемента, и силы, обусловленной начальной деформацией,
st
D
B
R
T
e
}
]{
[
]
[
}
{
0




(4.58)
Тетраэдральный симплекс-элемент. Выражение для матрицы упругости в трехмерном случае выглядит как




























0 0
0 0
0 0
0 0
0 0
0 0
0 0
0 0
0 0
1 0
0 0
1 0
0 0
1
]
[D
,
(4.59) где
)
2 1
)(
1
(
)
1
(
E








,
)
1
(





,
2
)
1
(
)
2 1
(






Такая запись подразумевает, что деформации объединены в вектор в следующем порядке
T
zx
yz
xy
z
y
x
T
}
{
}
{








Выражение для матрицы [B] получается дифференцированием соотно- шений (4.52) и состоит из четырех блоков
],
[
]
[
m
k
j
i
B
B
B
B
B

(4.60)

136 где














































x
N
z
N
y
N
z
N
x
N
y
N
z
N
y
N
x
N
B
i
i
i
i
i
i
i
i
i
i
0 0
0 0
0 0
0 0
0
]
[
(4.61)
Выражение для матрицы жесткости, определяемой соотношением, можно проинтегрировать точно, так как компоненты деформаций и напря- жений постоянны внутри элемента.
Подматрица матрицы жесткости имеет размерность 3х3 и определяется соотношением
e
s
T
r
e
rs
V
B
D
B
K
]
][
[
]
[
]
[

(4.62) где V
e
– объем тетраэдра; (r =1,…,4; s =1,…,4).
Узловые силы, обусловленные начальной деформацией, записываются в виде
e
T
e
V
D
B
R
}
]{
[
]
[
}
{
0




,
(4.63)
Трехмерный изопараметрический элемент. Особенностью определения матрицы жесткости данного элемента является разделение геометрических параметров элемента и физических параметров его материала. Формула для вычисления блока матрицы жесткости равна [4, 5]


V
s
T
r
e
rs
dV
B
D
B
K
]
][
[
]
[
]
[
,
(4.64) где матрица [B
r
] определяется соотношением (4.61), a [D] - выражением
(4.59); r =1,2,...,8.
Преобразуем (4.56) к следующему виду:

137
 








3 1
3 1
]
[
]
[
n
m
V
m
s
n
r
nm
e
rs
dV
x
N
x
N
D
K
(4.65)
Здесь x n
(n = 1,2,3) - глобальные координаты; [D
nm
] - матрицы размер- ности 3х3, полученные из матрицы [D] за счет коэффициентов, стоящих на пересечении n-х строк и m-х столбцов. Причем индекс n пробегает после- довательно три значения, соответствующие расположению
n
i
x
N


в каж- дой строке матрицы [B
r
]
T
, а m пробегает значения, соответствующие рас- положению m
i x
N


в столбцах матрицы [B
s
].
Интеграл вида






V
m
s
n
r
nm
rs
dV
x
N
x
N
K ]
[
,
(4.66) вычисляем, переходя к локальным координатам ξ, η, ζ
  
  





1 1
1 1
1 1
]
[



d
d
d
I
x
N
x
N
K
m
s
n
r
nm
rs
,
(4.67) где






































z
z
z
y
y
y
x
x
x
I ]
[
- матрица Якоби,
(4.68) где а[I]=det([I]) – якобиан матрицы.
Изопараметрический квадратичный элемент. При определении выра- жения для матрицы жесткости используется подход, описаний выше. В этом случае формулу для вычисления блока матрицы жесткости размерно- сти 2х2 можно представить в виде
 








2 1
2 1
]
[
]
[
n
m
V
m
s
n
r
nm
e
rs
dV
x
N
x
N
D
K
(4.69)

138
Здесь x
n
- декартовы координаты; [D
nm
] - матрицы размерности 2х2, полученные из матрицы [D].
При вычислении интегралов типа (4.67) в выражении (4.69) использу- ются локальные координаты ξ и η, а матрица Якоби имеет вид



















y
y
x
x
I ]
[
(4.70)
Одно из главных преимуществ формул (4.66) и (4.56) проявляется при решении физически нелинейных задач, когда изменяется только матрица упругости. Интегралы типа (4.67) вычисляются численно с использовани- ем квадратур Гаусса [5]:
- для плоского элемента
 
 

1 1
1 1
)]
,
(
[




d
d
G
I
,
(4.71)
- для объемного элемента
  
  

1 1
1 1
1 1
)]
,
,
(
[






d
d
d
G
I
,
(4.72) где G – некоторая функция, выражения для которой зависят от типа эле- мента, числа точек интегрирования и даны в [5].
Отметим, что, используя формулы (4.66) и (4.69), не представляет тру- да рассчитать матрицы жесткости и для элементов более высокого поряд- ка.
Имея матрицы жесткости отдельных элементов, можно получить гло- бальную матрицу жесткости рассматриваемой области [K
e
]. Для получения используется процедура объединения по элементам [5,7], то есть элемент- ная матрица [К]
е
прибавляется к системной матрице [К
е
] сразу же после вычисления. Вектор-столбец узловых сил получается аналогичным обра- зом.

139
Обычное решение
i
i
i
i
i
i
i
i
i
d
d
d
}
{
}
{
}
{
}
{
}
{
}
{
}
{
}
{
}
{
1 1
1


















, полученное МКЭ, является приближением к истинному или точному решению теории упругости.
Сформулируем критерии сходимости, которые непосредственно сле- дуют из сказанного выше, и состоят в следующем:
1. Для того чтобы выполнялось требование сходимости, необходимо, чтобы представления перемещения и любой его производной, появляю- щейся в функционале, стремились к их точным значениям на каждом эле- менте по мере стремления к нулю размера элемента.
2. Условие сходимости состоит в том, что по мере стремления к нулю размера элемента члены с производными и функцией в функционале должны стремиться к функции той же гладкости, что и точное решение
(предполагается непрерывность точного решения).
Критерии 1 (полноты) и 2 (согласованности) являются достаточными условиями сходимости вариационного метода конечных элементов. При условии выполнения критериев полноты и согласованности путем доста- точного уменьшения размера элемента теоретически можно достичь лю- бой требуемой точности решения.
4.4 Другие формулировки МКЭ
Несмотря на то, что приближенная минимизация функционала - самый распространенный способ подхода к МКЭ, это никоим образом не означа- ет, что такой подход является единственно возможным.
Наиболее популярными из других вариантов МКЭ являются методы, позволяющие получить основные соотношения МКЭ непосредственно из дифференциальных уравнений задачи; возможные преимущества таких методов состоят в том, что: а) исчезает необходимость искать функциональный эквивалент извест- ным ДУЧП; б) эти методы могут быть распространены на задачи, для которых функционал - либо вообще не существует, либо пока еще не получен.

140
Рассмотрим задачи приближенного решения системы ДУЧП, которой должна удовлетворять неизвестная функция в области V . Запишем основ- ное уравнение в виде
0
)
(


L
,
(4.73) а граничное условие на границе С как
0
)
(


C
(4.74)
Если пробная функция, удовлетворяющая граничным условиям, запи- сана в общей форме (4.3)
 
 


N

ˆ
,
(4.75) где [N]- матрица базисных функций; {

} - основные неизвестные задачи, то в общем случае
 
0
ˆ



R
L

,
(4.76) где
)
ˆ
(
)
(


L
L
R


- невязка решения.
Наилучшим решением будет то, которое дает во всех точках области V наименьшую невязку R .
Очевидно, что решение можно получить, использовав то обстоятель- ство, что если невязка R тождественно равна нулю всюду в области, то
0


v
WRdV
,
(4.77) где W – любая функция координат. Если число неизвестных параметров
{

} равно n, то, выбрав n линейно независимых функций W
i
, запишем со- ответствующую систему уравнений

141
 
 
0
)
(




dV
N
L
W
WRdV
v
i
v

,
(4.78) из которой может быть найдена сеточная функция {

}.
Описанный процесс называется методом взвешенных невязок, a W
i
- весовой функцией. Выбор различных W
i приводит к различным классиче- ским вариантам МКЭ. Наиболее популярным из них является метод Галер- кина. В этом случае W
I
=N
I
, т.е. в качестве весовой функции выбирается функция формы, с помощью которой аппроксимируется решение {

}. Этот метод обычно приводит к наилучшим результатам.
Отметим один недостаток метода взвешанных невязок. В этом методе дифференциальный оператор L, содержит производные более высоких по- рядков, чем вариационный функционал F. Поэтому необходимо обеспе- чить выполнение условий непрерывности функций формы более высокого порядка.
4.5 Классификация конечных элементов
Наиболее очевидная классификация элементов делит их:
- на одномерные;
- двумерные;
- трехмерные.
Эти группы могут разделяться в зависимости от того, включают ли уз- ловые переменные только значение функций (лагранжевые элементы) или также и значение производных (эрмитовы элементы).
4.6 Базисные функции конечных элементов и выбор конечного элемента
Произвольная пробная функция, определенная на элементе, записыва- ется в линейной форме
 
 
e
e
e
U
N
U

ˆ
,
(4.79) где [N]-матрица базисных функций;
{U}
e
-узловой вектор элемента.

142
Для лагранжева элемента в узле имеется только одна степень свободы- значение функции, и, следовательно, уравнение (4.79) может быть записа- но в виде























s
k
s
k
U
U
U
U
N
N
N
N
U


2 1
2 1
ˆ
,
(4.80) где 1,…, K ,…, S-номера узлов, а S-общее количество узлов элемента.
Для эрмитова элемента каждая из компонент U
k должна записываться как столбец, включающий и производные указанной функции, которые в этом случае также являются узловыми переменными. Если каждый из S узлов элемента e имеет q степеней свободы, то каждая компонента U
k в уравнении (4.80) представляет собой столбец
 















kq
k
k
k
U
U
U
U

2 1
,
(4.81) где k=1,2,…,S и аналогичные базисные функции [N
k
] должны записываться как строки
 


kq
k
k
k
N
N
N
N

2 1

,
(4.82)
Таким образом, уравнение (4.80) при соответствующем выборе переменных является общим для случаев лагранжевых эрмитовых эле- ментов.
Ранее мы с вами установили, что базисные функции являются функци- ями независимых переменных (X,Y) и узловых координат. Для получения базисных функций любого отдельного элемента применимы два подхода:
-в одном из них используются обобщенные координаты;

143
-в другом интерполяционные формулы.
4.6.1 Получение базисных функций в обобщенных координатах
Этот подход особенно удобен для простых элементов, использующих полные полиномы низкого порядка. В случае более сложных элементов указанный подход становится неэффективным, и вместо него обычно при- меняют метод интерполяции.
Рассмотренный подход ранее был проиллюстрирован на примере тре- угольного-симплекс-элемента. Применим его для прямоугольного элемен- та со сторонами параллельными осям глобальной системы координат OXY
(рисунок 4.7).
Простейшая пробная функция для элемента содержит только четыре неизвестных параметраL
i