1329
.pdf= 1
4µSтр
|
1 ∂ |
|
∂ |
A 2 |
|
∂ |
A 2 |
|
|
|
||||||||||
|
|
|
|
|
|
∫∫ |
|
|
|
+ |
|
|
|
dxdy = |
|
|
||||
|
2µ ∂ |
A |
|
|
|
y |
|
|
||||||||||||
|
|
|
∂ |
x |
|
∂ |
|
|
|
|
|
|||||||||
|
|
|
|
|
n |
|
|
|
|
|
|
|
|
|
|
(6.79) |
||||
|
|
|
|
|
|
S тр |
|
|
|
|
|
|
|
|
|
|
|
|
||
|
|
+ |
Ambm |
+ |
Anbn )bn |
+ |
( Alcl |
+ |
Amcm |
+ |
. |
|||||||||
( Albl |
|
|
|
|
|
|
|
|
Ancn )cn |
Если принять, что в каждом элементе плотность стороннего тока является постоянной величиной, то
|
|
|
|
|
∂ |
( J |
|
A) = J |
|
∂ A |
; |
|
|
|
∂ |
( J |
|
A) = J |
|
∂ A |
; |
|
|||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||||||||||
|
|
|
|
|
∂ A |
ст |
|
|
ст∂ A |
|
|
|
|
∂ A |
|
|
|
|
ст |
|
ст∂ A |
|
|||||||
|
|
|
|
|
l |
|
|
|
|
|
|
l |
|
|
|
|
m |
|
|
|
|
m |
|
||||||
|
|
|
|
|
|
|
|
|
∂ |
|
( J |
|
A) = J |
|
|
∂ |
|
A . |
|
|
|
(6.80) |
|||||||
|
|
|
|
|
|
|
|
|
|
|
|
ст∂ |
|
|
|
|
|||||||||||||
|
|
|
|
|
|
|
|
|
∂ A |
ст |
|
|
|
|
|
|
A |
|
|
|
|
|
|||||||
|
|
|
|
|
|
|
|
|
|
n |
|
|
|
|
|
|
|
|
|
|
|
n |
|
|
|
|
|
||
Тогда |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||
|
∂ |
|
|
∫∫ Jст Adxdy = Jст |
∫∫ Adxdy = Jст |
1 |
|
∫∫ (al + blx + cl y)dxdy ; |
(6.81) |
||||||||||||||||||||
|
∂ A |
|
2S |
||||||||||||||||||||||||||
|
l S тр |
S тр |
|
|
|
|
|
|
|
|
тр |
S тр |
|
|
|
|
|
|
|
||||||||||
|
∂ |
|
∫∫ Jст Adxdy = Jст ∫∫ Adxdy = Jст |
|
1 |
∫∫ (am + bmx + cm y)dxdy ; |
(6.82) |
||||||||||||||||||||||
|
∂ A |
2S |
|||||||||||||||||||||||||||
|
m S тр |
S тр |
|
|
|
|
|
|
|
|
тр |
S тр |
|
|
|
|
|
|
|
||||||||||
|
∂ |
|
|
|
|
|
|
|
|
|
|
|
1 |
|
|
∫∫ (an + bnx + cn y)dxdy . |
|
||||||||||||
|
|
∫∫ Jст Adxdy = Jст ∫∫ Adxdy = Jст |
|
(6.83) |
|||||||||||||||||||||||||
|
∂ An |
2Sтр |
|||||||||||||||||||||||||||
|
|
|
|
S тр |
S тр |
|
|
|
|
|
|
|
|
|
|
|
|
S тр |
|
|
|
|
|
В записанных выражениях функции под знаком интеграла зависят от пространственных координат. Поэтому вычисление интеграла необходимо производить по площади треугольника. Положим для простоты вычислений, что рассматривается прямоугольный треугольник, изображённый на рис. 6.2.
151
Рассмотрим вычисление |
интеграла |
|
в выражении (6.83). В данном случае |
||
координата х изменяется в |
пределах |
|
xl ≤ x≤ |
xm . По координате y рассматри- |
|
ваемый |
треугольник ограничен прямой |
|
y = yl |
снизу и прямой y = yl |
+ k (x − xl ) |
сверху. |
Поэтому |
|
Рис. 6.2. К определению |
|
(6.84) |
интеграла (6.83) |
|
|
Возьмём интеграл
(6.85)
an = xl ym − xm yl = xl ym − xm yl + xl yl − xl yl =
= xl ( ym − yl ) + yl (xl − xm) = xl 0 + yl (−hx) = − yl hx;
bn = yl − ym = 0; |
cn = xm − xl = h x; |
k = |
y n − yl |
= |
h y |
. |
||
xn |
− xl |
h x |
||||||
|
|
|
|
|
Тогда
(6.86)
(6.87)
(6.88)
152
x |
( x − xl ) |
2 |
dx = h x k 2 ( x − xl ) |
3 |
|
|
|
||||
∫m h x k 2 |
|
|
|||
xl |
2 |
|
2 |
3 |
|
|
|
|
|
|
xm
= h x k 2h3x . |
(6.89) |
6 |
|
xl
Подставляя в полученный результат k согласно (6.87) и учитывая, что площадь прямоугольного треугольника
S тр = 1 h x h y , 2
в окончательном виде будем иметь
∂ |
∫∫ J ст Adxdy = |
|
1 |
h x h2y h3x |
J ст = |
S трJ ст |
. |
(6.90) |
|
∂ A |
2 |
S тр |
6 |
2 |
3 |
||||
m |
S тр |
|
|
h x |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Остальные интегралы (6.81) и (6.82) вычисляются аналогично и имеют те же значения.
Таким образом, минимум функционала (6.68) может быть записан в виде следующей системы уравнений:
1 |
|
bl2 |
+ cl2 |
|
|
||||
|
blbm |
+ cl cm |
||
|
|
|||
4µS тр |
||||
blbn + cl cn |
blbm + cl cm |
blbn + cl cn |
|
Al |
|
|
|
1 |
|
|
|
|
|
|
||||||
bm2 + cm2 |
bmbn + cmcn |
× |
Am |
= |
S трJ ст |
|
1 |
. (6.91) |
|
3 |
|||||||||
bmbn + cmcn |
bn2 + cn2 |
|
An |
|
|
1 |
|
||
|
|
|
|
Указанная система уравнений записывается для всех элементов исследуемой области, а её решение с учётом известных значений векторного потенциала на границе области позволяет рассчитать искомые значения векторного потенциала во всех узлах исследуемой области.
Указанная система алгебраических уравнений характеризует поэлементное объединение, процедура которого сводится к вычислению матриц для каждого элемента.
Помимо поэлементного объединения возможно объединение по узлам, когда уравнение записывается для каждого узла исследуемой
153
области. Система алгебраических уравнений в этом случае может быть получена исходя из следующих соображений.
Каждое из уравнений системы (6.91) получено путём минимизации энергетического функционала (6.68). При этом для получения минимума производилось дифференцирование функционала по значениям векторного потенциала в узлах треугольника, и полученные производные приравнивались к нулю. В исследуемой области каждый узел может принадлежать одновременно нескольким элементам. Система алгебраических уравнений, записанная для каждого элемента этой группы, будет содержать, естественно, производные по значению векторного потенциала рассматриваемого узла. В этом случае можно, просуммировав уравнения рассматриваемого узла для всех элементов, включающих в себя данный узел, получим для него уравнение, которое содержит значения векторного потенциала в соседних узлах.
Для примера рассмотрим простейшую краевую задачу для уравнения Пуассона в прямоугольной области с нулевыми граничными условиями. Разобьём рассматриваемую область на прямоугольные ячейки со сторонами величиной h x, h y и вершинами, совпадающими
с узлами сетки. Каждую ячейку области разобьём диагональю на два прямоугольных треугольника с катета-
ми, равными hx и h y . В результате об-
ласть оказывается разбитой на конечное число прямоугольных треугольни-
ков – конечных элементов.
Рассмотрим узел №1 этой области, являющейся вершиной элементарных треугольников, обозначенных цифрами
1–2–3–4–5–6 ( рис. 6.3).
Векторный потенциал для вершин каждого из прямоугольных треугольников описывается системой уравнений (6.91). Для узла № 1 имеем
Рис. 6.3. Квыводууравнения векторного потенциала в узле № 1
154
|
∂ F |
|
6 |
|
F |
|
|
|
|
|
= ∑ ∂ |
, |
(6.92) |
||||
|
|
|||||||
|
∂ A1 ∑ |
i=1 |
∂ |
A1 i |
|
где i − номер треугольника с вершиной в узле № 1; F – функционал, определяемый выражением (6.68).
Производная по А1 в каждом элементе соответствует первой
∂ F
строке системы (6.91). Определим для первого треугольника
∂ A1 1
(см. рис. 6.3). Обозначим его вершины, принимая во внимание, что обход производится против часовой стрелки: l = 1; m = 3; n = 2. Тогда в соответствии с (6.62), (6.63)
bl = y3 − y2 = 0; |
bm = y2 − y1 = −h y; |
bn = y1 − y3 = h y; |
(6.93) |
cl = x2 − x3 = h x; |
cm = x1 − x2 = 0; cn = x3 − x1 = −h x. |
(6.94) |
|
Приведём уравнения системы (6.91) |
к общему знаменателю |
и определим правые части уравнений:
4µS тр |
S трJ ст |
|
J ст |
|
1 |
2 2 |
J ст |
2 2 |
|
|
= 4µ |
|
|
hx h y = µ |
|
h xh y . |
|||
3 |
3 |
4 |
3 |
||||||
|
|
|
|
|
Первое уравнение из этих уравнений
|
∂ F |
|
= (bl2 + cl2 ) A1 + (blbm + clcm ) A3 + |
|
|
|
|||
∂ A1 |
||||
|
1 |
|
+(blbn + cl cn ) A2 = µ J ст h2xh2y. 3
(6.95)
(6.96)
Подставляя в это уравнение значения bl , bm , bn , cl , cm , cn и выполняя преобразования, получим
155
|
∂ F |
|
= h2x A1 − h2x A2 |
= µ |
J ст |
|
|
|
|
|
h2xh2y . |
(6.97) |
|||||
∂ A1 |
|
|||||||
|
1 |
|
3 |
|
|
Аналогичные операции проделаем для всех треугольников с вер-
шиной в узле № 1. Полученные данные сведём в табл. 6.1. |
|
|
|
||||||||
|
|
|
|
|
|
|
Таблица |
6 . 1 |
|||
Коэффициенты конечных элементов для узла № 1 |
|
||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
Номер |
l |
m |
n |
bl |
bm |
bn |
cl |
|
cm |
|
cn |
треугольника |
|
|
|||||||||
|
|
|
|
|
|
|
|
|
|
|
|
1 |
1 |
3 |
2 |
0 |
–hy |
hy |
hx |
|
0 |
|
– hx |
2 |
1 |
4 |
3 |
hy |
– hy |
0 |
0 |
|
hx |
|
– hx |
3 |
1 |
5 |
4 |
hy |
0 |
– hy |
– hx |
|
hx |
|
0 |
4 |
1 |
6 |
5 |
0 |
hy |
– hy |
– hx |
|
0 |
|
hx |
5 |
1 |
7 |
6 |
– hy |
hy |
0 |
0 |
|
– hx |
|
hx |
6 |
1 |
2 |
7 |
– hy |
0 |
hy |
hx |
|
– hx |
|
0 |
Врезультате получим длявсех элементов следующие уравнения:
|
|
∂ F |
|
= h2y A1 − h2y A4 |
= µ |
|
J ст |
|
|
|
|
|
|
|
|
||||
|
|
|
|
|
h2xh2y ; |
|
|
|
|
|
(6.98) |
||||||||
|
∂ A1 |
|
|
|
|
|
|
|
|
||||||||||
|
|
2 |
|
3 |
|
|
|
|
|
|
|
|
|
||||||
|
∂ F |
|
|
= (h2x + h2y) A1 − h2x A5 − h2y A4 |
= µ |
J ст |
|
|
|
||||||||||
|
|
|
h2x h2y ; |
(6.99) |
|||||||||||||||
∂ A1 |
|
|
|||||||||||||||||
|
3 |
|
|
|
|
|
|
|
|
3 |
|
|
|
||||||
|
∂ F |
|
|
= h2x A1 − h2x A5 |
= µ |
|
J ст |
|
|
|
|
|
|
|
|
|
|||
|
|
|
|
h2xh2y ; |
|
|
|
|
|
(6.100) |
|||||||||
∂ A1 |
|
|
|
|
|
|
|
||||||||||||
|
4 |
|
3 |
|
|
|
|
|
|
|
|
|
|||||||
|
∂ F |
|
|
= h2y A1 − h2y A7 |
= µ |
J ст |
|
|
|
|
|
|
|
|
|||||
|
|
|
h2x h2y ; |
|
|
|
|
|
(6.101) |
||||||||||
∂ A1 |
|
|
|
|
|
|
|
||||||||||||
|
5 |
|
3 |
|
|
|
|
|
|
|
|
|
|||||||
|
∂ F |
|
|
= (h2x + h2y ) A1 − h2x A2 − h2y A7 |
= µ |
J ст |
|
|
|||||||||||
|
|
h2xh2y . |
(6.102) |
||||||||||||||||
∂ A1 |
|
||||||||||||||||||
|
6 |
|
|
|
|
|
|
|
|
3 |
|
|
|
156
Суммируя левые и правые части этих уравнений, получим
|
∂ F |
|
= 4(hx2 + hy2 ) A1 − 2hx2 A2 |
− 2hy2 A4 |
− 2hx2 A5 |
− |
|
|
|
||||||
∂ A1 |
|||||||
|
Σ |
|
|
|
(6.103) |
− 2hy2 A7 − 2µJстhx2 hy2 .
Учитывая равенство нулю производных в левой части уравнения, выполняем элементарные преобразования и получаем уравнение
A4 − 2 A1 + A7 |
+ |
A2 − 2 A1 + A5 |
= −µJ ст , |
(6.104) |
h2x |
|
|||
|
h2y |
|
совпадающее суравнением, полученным конечно-разностным методом. В качестве примера рассмотрим уравнение Лапласа
∂ 2u |
+∂ |
|
2u |
= 0 |
, |
(6.105) |
∂ x2 |
|
|
||||
∂ |
y2 |
|
|
решаемое методом конечных элементов в работе [35] в прямоугольной области [0:2, 0:2] с краевыми условиями:
u(x,0) = 50; u(x, 2) =100 ; |
∂ u |
(0, y) = 0 ; |
∂ u |
(2, y) = 0 . |
|
|
|||
|
∂ x |
∂ x |
Рассматриваемая область разбита на 16 треугольников с 15 узлами (рис. 6.4).
Система алгебраических уравнений, соответствующая краевой задаче (6.105), получена в работе [35], записана в матричном виде (6.107), а её решение представлено ниже:
u1 = 50, u2 |
= 50, u3 = 50 − краевые условия. |
u4 = 62,5, |
u5 = 62,5, u6 = 62,5, u7 = 75, u8 = 75, |
157
u9 |
= 75, u10 = 87,5, |
u11 = 87,5, u12 = 87,5, |
u13 |
=100, u14 =100, |
u15 =100 − краевые условия. |
Рис. 6.4. Триангуляция исследуемой области [0:2, 0:2]
Покажем, что каждое уравнение из рассматриваемой системы соответствует конечно-разностной аппроксимации дифференциальных операторов решаемого уравнения Лапласа.
Рассмотрим, например, уравнение для узла № 8, предварительно разделив его коэффициенты и правую часть на 2:
−4u5 −1u7 +10u8 −1u9 − 4u11 = 0 . |
(6.106) |
158
1 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
|
|
|
0 |
0 |
0 |
|
0 |
0 |
|
0 |
|
u1 |
|
50 |
|
||
0 |
1 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
|
|
|
0 |
0 |
0 |
|
0 |
0 |
|
0 |
|
u2 |
|
50 |
|
||
0 |
0 |
1 |
0 |
0 |
0 |
0 |
0 |
0 |
|
|
|
0 |
0 |
0 |
|
0 |
0 |
|
0 |
|
u3 |
|
50 |
|
||
−4 0 |
0 |
10 −2 0 −4 0 |
0 |
|
|
0 0 0 0 0 0 |
|
u4 |
|
0 |
|
|||||||||||||||
0 |
−8 0 −2 20 −2 |
0 |
−8 0 |
|
0 |
0 |
0 |
|
0 |
0 |
|
0 |
|
u5 |
|
0 |
|
|||||||||
0 |
0 |
−4 0 −2 10 0 0 −4 0 |
0 |
0 |
|
0 |
0 |
|
0 |
|
u6 |
|
0 |
(6.107) |
||||||||||||
0 |
0 |
0 |
−4 0 |
0 |
10 −2 0 −4 0 |
0 |
0 0 0 |
|
u7 |
|
0 |
|||||||||||||||
|
|
|
||||||||||||||||||||||||
0 |
0 |
0 |
0 |
−8 0 −2 20 −2 |
|
0 |
−8 0 |
|
0 |
0 |
|
0 |
|
u8 |
= |
0 |
|
|||||||||
0 |
0 |
0 |
0 |
0 |
−4 0 −2 10 0 0 −4 0 |
0 |
|
0 |
|
u9 |
|
0 |
|
|||||||||||||
0 |
0 |
0 |
0 |
0 |
0 |
−4 0 |
0 |
|
|
10 −2 0 −4 0 |
0 |
|
u10 |
|
0 |
|
||||||||||
0 |
0 |
0 |
0 |
0 |
0 |
0 |
−8 0 −2 20 −2 |
|
0 |
−8 0 |
|
u11 |
|
0 |
|
|||||||||||
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
−4 |
|
0 |
−2 10 0 0 −4 |
|
u12 |
|
0 |
|
||||||||||
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
|
|
|
0 |
0 |
0 |
|
1 |
0 |
|
0 |
|
u13 |
|
100 |
|
||
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
|
|
|
0 |
0 |
0 |
|
0 |
1 |
|
0 |
|
u14 |
|
100 |
|
||
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
|
|
|
0 |
0 |
0 |
|
0 |
0 |
|
1 |
|
u15 |
|
100 |
|
||
|
Уравнение Лапласа в конечно-разностном виде |
|
|
|
|
|
|
|||||||||||||||||||
|
|
|
|
ui +1, j − 2ui, j + ui −1, j |
+ |
ui, j +1 − 2ui, j + ui, j −1 |
= 0 . |
|
|
|
||||||||||||||||
|
|
|
|
|
|
|
|
|
|
|
||||||||||||||||
|
|
|
|
|
|
|
h2x |
|
|
|
|
|
|
|
|
h2y |
|
|
|
|
|
|
|
|
|
|
|
Тогда для узла № 8 при hx =1,0 и h y = 0,5 получим |
|
|
|
||||||||||||||||||||||
|
|
|
|
|
|
u9 − 2u8 + u7 |
+ |
u11 − 2u8 + u5 |
= 0 . |
|
|
|
|
|
|
|
||||||||||
|
|
|
|
|
|
|
|
2 |
|
|
|
|
|
0,5 |
2 |
|
|
|
|
|
|
|
|
|
||
|
|
|
|
|
|
|
|
1 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Выполнив элементарные преобразования, получим уравнение
(6.106).
В ряде случаев энергетический функционал не известен или его нахождение связано со значительными трудностями. В этом случае система алгебраических уравнений может быть получена непосредственно при интегрировании уравнений Максвелла. Рассмотрим процедуру получения системы алгебраических уравнений для этого случая.
159
Согласно определению ротор любого вектора может быть выражен через циркуляцию следующим образом:
(6.108)
Введём векторный потенциал rot A = B и запишем уравнение Максвелла в виде
rot |
1 |
rot |
|
= |
|
. |
(6.109) |
|
A |
J |
|||||||
µ |
||||||||
|
|
|
|
|
|
|
||
Интегрируя полученное |
выражение по произвольной |
площади |
ииспользуя определение ротора, этому выражению можно придатьвид
(6.110)
Таким образом, решение дифференциального уравнения (6.109) можно свести к системе уравнений (6.110), записываемых для всех участков площади исследуемой области.
Положим, что исследуемая область разбита на определённое число элементарных площадок в виде треугольников (триангулирована), в каждом из которых векторный потенциал является линейной функцией пространственных координат и записывается в виде
A = A1 |
− |
x − x1 |
|
y2 − y1 |
A2 − A1 |
|
+ |
y − y1 |
|
|
x2 − x1 |
A2 |
− A1 |
|
= |
|
|
|
|
|
|||||||||||||
S |
|
y3 − y1 |
A3 − A1 |
|
S |
|
x3 |
− x1 |
A3 − A1 |
|
||||||
|
|
|
|
|
|
|
|
(6.111)
= A1 − K1 ( x − x1 ) + K2 ( y − y1 ),
где коэффициенты
K 1 |
= |
1 |
|
|
y |
2 − y1 |
A2 − A1 |
|
; |
K 2 |
= |
1 |
|
x2 − x1 |
A2 |
− A1 |
|
. (6.112) |
|
|
|
|
|
||||||||||||||||
S |
|
y3 |
− y1 |
A3 − A1 |
|
S |
x3 − x1 |
A3 − A1 |
|
||||||||||
|
|
|
|
|
|
|
|
|
160