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

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 соответственно, то дифференцирова-