Подставляя теперь (28), (34) и (35) в формулы (20), (21), получаем
разрешающую систему для КЭ в матричной форме:
[
K
(
q
)]
12
×
12
{
q
}
12
− {
N
p
(
q
)
}
12
×
3
{
y
}
3
=
{
P
}
12
;
[
B
p
(
q
)]
т
3
×
12
{
q
}
12
=
{
h
p
(
q
)
}
3
,
(36)
где обозначено:
[
K
(
q
)]
12
×
12
=
R
◦
V
([
B
]
т
+ [
a
B
]
т
+ [ ˜
B
]
т
)[
A
]([
B
] + [
a
B
])
d
◦
V
;
[
B
p
(
q
)]
12
×
3
=
Z
◦
V
([
B
]
т
+ [
a
B
]
т
)[Φ
p
]
d
◦
V
;
[
N
p
(
q
)]
12
×
3
=
Z
◦
V
([
B
]
т
+ [
a
B
]
т
+ [ ˜
B
]
т
)([Φ
p
] + [
I
]([
B
] + [
a
B
])[Φ
q
]
{
y
}
)
d
◦
V
;
{
P
}
=
Z
◦
Σ
[Φ]
т
{
S
}
d
◦
S
;
{
h
p
(
q
)
}
=
Z
◦
V
h
{
Φ
p
}
d
◦
V .
(37)
Здесь также обозначены матрицы, полученные тензорным умножени-
ем столца на строку:
[Φ
p
]
4
×
3
=
{
I
0
}
4
{
Φ
p
}
т
3
; [Φ
q
]
15
×
3
=
{
q
}
15
{
Φ
p
}
т
3
.
(38)
Численный алгоритм.
На основе локальной системы алгебраиче-
ских уравнений (37), записанной для одного КЭ, далее была построе-
на глобальная система для всей области
◦
V
. Обе эти системы являются
нелинейными, поскольку матрица
[
a
B
]
(см. (27)), а следовательно, и
матрицы
[
K
]
,
[
B
p
]
,
[
N
p
]
зависят от вектора неизвестных
{
q
}
.
Для решения системы нелинейных уравнений применяли итераци-
онный метод, при котором на каждом текущем шаге итерации значения
компонент вектора неизвестных
{
q
}
в матрицах
[
K
]
,
[
B
p
]
,
[
N
p
]
и век-
торах
{
P
}
и
{
h
p
}
выбирали с предыдущего шага итерации. Решение
линеаризованных систем было проведено методом Гаусса.
Генерацию конечно-элементной (КЭ) треугольной сетки осуще-
ствляли “вручную”, с помощью преобразования прямоугольника в
двумерную область решения задачи
◦
V
в координатах
r, z
, при этом
полагали, что
◦
V
ограничена двумя торцевыми плоскостями
z
= 0
и
z
=
z
max
(см. далее рис. 4,
а
). В результате была получена симметрич-
ная относительно плоскости
z
=
1
2
z
max
КЭ сетка, которая, как показали
76
ISSN 1812-3368. Вестник МГТУ им. Н.Э. Баумана. Сер. “Естественные науки”. 2007. № 3
