Электромагнитные методы
УДК 620.179.14
РАСЧЕТ НАПРЯЖЕННОСТИ МАГНИТНОГО ПОЛЯ ОТ ОДНОРОДНОГО
ЦИЛИНДРА КОНЕЧНЫХ РАЗМЕРОВ, ПОМЕЩЕННОГО В ПРОИЗВОЛЬНОЕ
ВНЕШНЕЕ ПОЛЕ
© 2022 г. В.В. Дякин1, О.В. Кудряшова1,*, В.Я. Раевский1,**
1Институт физики металлов имени М.Н. Михеева УрО РАН,
Россия 620137 Екатеринбург, ул. С. Ковалевской, 18
E-mail: *kudryashova_ov@imp.uran.ru; **ravskii@mail.ru
Поступила в редакцию 09.03.2022; после доработки 20.03.2022
Принята к публикации 22.03.2022
Представлен алгоритм нахождения напряженности результирующего магнитного поля внутри и вне однородного
цилиндра конечных размеров, помещенного во внешнее магнитное поле произвольной конфигурации, сводящийся в
основной его части к решению некоторого количества систем трех одномерных линейных интегральных уравнений.
Описаны сложности и способы их преодоления при дискретизации этих уравнений для соответствующей программной
реализации. Проведено сравнение результатов работы реализующей описанный алгоритм компьютерной программы в
предельных частных случаях с известными аналитическими ответами.
Ключевые слова: интегродифференциальное уравнение магнитостатики, цилиндрический магнетик, магнитный не-
разрушающий контроль.
DOI: 10.31857/S0130308222040078, EDN: BLOMAU
1. ВВЕДЕНИЕ
Для решения многих практических задач из области магнетизма, например, задач неразруша-
ющего магнитного контроля, актуальной является проблема создания обоснованных алгоритмов
аналитического или численного решения задач магнитостатики по вычислению напряженности
результирующего поля применительно к магнитным телам различной формы, помещенным во
внешнее магнитное поле. Обоснованные эффективные методы решения такого рода задач имеются
по большей части для безграничных модельных осесимметричных тел относительно простой гео-
метрической формы, помещенных в постоянное внешнее магнитное поле определенного (удобного
для аналитического исследования) направления, что позволяло во многих случаях свести задачу к
двумерной и пренебречь краевыми эффектами. Однако большой теоретический и практический
интерес представляют задачи для реальных тел конечных размеров в произвольном внешнем поле.
Что касается тел с цилиндрической симметрией, давно исследована задача для бесконечно
длинного цилиндра с постоянной магнитной проницаемостью, помещенного в бесконечную маг-
нитную среду с иной магнитной проницаемостью, в однородном внешнем поле, решение которой
представляется в элементарных функциях [1, с. 245]. В [2] изучена задача расчета напряженности
результирующего поля бесконечного магнитного цилиндра при условии неоднородного намагни-
чивания, решение которой записывается через специальные функции. В [3] эта задача исследована
для однородного цилиндра конечных размеров в произвольном внешнем поле. В этой работе к
возникающему двумерному интегральному уравнению (интегрирование по полной поверхности
цилиндра) непосредственно применен метод коллокаций, что приводит к системе линейных урав-
нений достаточно большой размерности, а для более комфортного вычисления матричных элемен-
тов принят ряд упрощающих предположений.
В статье авторов [4] предложен и обоснован подход, сводящий вычисление напряженности маг-
нитного поля от однородного цилиндра конечных размеров, помещенного в произвольное внешнее
поле, к решению некоторого количества систем трех одномерных линейных интегральных урав-
нений. Настоящая работа является продолжением этой статьи и посвящена тонкостям и результа-
там программной реализации этого подхода, тестированию получаемых результатов на их соот-
ветствие физическим законам явления, а также всесторонней проверке адекватности получаемых
результатов на предельных случаях с известными «аналитическими» ответами.
Создание подобных алгоритмов и программ расчетов полей с контролируемой точностью дик-
туется также необходимостью тестирования известных пакетов универсальных программ (типа
ELCUT, ANSYS, ELMER), формальное использование которых приводит ко многим проблемам.
Подробному описанию недостатков этих программ и подводных камней при их использовании по-
священа работа [5].
64
В.В. Дякин, О.В. Кудряшова, В.Я. Раевский
2. СХЕМА РАСЧЕТА НАПРЯЖЕННОСТИ ПОЛЯ РЕАКЦИИ
Для решения задач магнитостатики мы исходим из так называемого основного уравнения маг-
нитостатики, которое в случае однородного магнетика с постоянной магнитной проницаемостью
µ имеет вид:
µ-1
H(r′)
0
H r)−
∇div
dr′
=
H r),
r∈R3\S.
(1)
∫
4π
| r
−r′|
Ω
Это уравнение эквивалентно системе уравнений Максвелла для случая магнитостатики
(см. [6, с. 17], [7, с. 149], [8]) и связывает искомую напряженность результирующего магнитного
поля H(r) ={Hx(r), Hy(r), Hz(r)} в произвольной точке пространства r =(x, y, z) (не лежащей на
0
0
0
0
границе магнетика) с напряженностью
H r)
={H r),
H r),H
(r)}
заданного поля внешнего
x
y
z
источника. В уравнении (1) Ω есть область в пространстве R3, ограниченная поверхностью S и за-
нятая исследуемым магнетиком с заданной постоянной магнитной проницаемостью µ.
Рассмотрим магнетик в форме прямого кругового цилиндра (см. рис. 1). Область Ω , занятая
магнетиком, представляет собой прямой круговой цилиндр длины l с боковой поверхностью S1, ось
которого совмещена с осью z, а нижнее и верхнее основания S2 и S3 суть круги радиуса R, располо-
женные в плоскостях z = d и z = d + l соответственно.
z
n
H0
S3
d + l
μ
S1
n
R
S2
d
n
y
x
Рис. 1. Однородный цилиндр в магнитном поле.
В нашей работе [4] разработана и обоснована методика вычисления так называемого поля реак-
ции (магнетика) HR(r):= H(r) - H0(r), которая сводится к схеме, представленной ниже.
R
R
R
R
Цилиндрические координаты поля реакции
H ϕ,z)
={
H r,ϕ,z),
H r,ϕ,z),
H r,ϕ,z)
}
в
r
ϕ
z
произвольной точке пространства с цилиндрическими координатами (r, φ, z), не лежащей на по-
верхности магнетика
S=S
S
S
, вычисляются по формулам:
1
2
3
Дефектоскопия
№ 4
2022
Расчет напряженности магнитного поля от однородного цилиндра конечных размеров...
65
∞
1
1
R
(1)
(2)
r
H r,ϕ,z)
= π(µ-1)
∑
cosnϕ[
t∫a
n
(z′
n
)α r, 1,
z -tz′)dz′ + ∫a r′
n
n
α r,
r′,
z)
dr′ +
n=0
0
0
1
1
(3)
(1)
+∫a r′
n
n
α r,
r′,
z-t)
dr′ +sinnϕ[
t∫b
n
(z′
n
)α r, 1,
z-t⋅z′)
dz′ +
0
0
1
1
(2)
(3)
+∫b r′
α r,
r′,
z
)
dr′ + ∫b r′
α
(r,
r′, z − t) dr′];
(2)
n
n
n
n
0
0
∞
1
R
(1)
H r,ϕ,z)
=- π(µ-
1)
cosnϕ[
t∫b
(z′)β
(r, 1,
z -tz′)dz
′+
ϕ
∑
n
n
n=1
0
1
1
1
(2)
(3)
(1)
+∫b r′
β r,
r′,
z)
dr′ + ∫b r′
β r,
r′,
z-t)
dr′] sinnϕ[t∫a
(z′)β r, 1,
z-t⋅z′)
dz′ +
n
n
n
n
n
n
0
0
0
1
1
(2)
(3)
+∫a r′
β r,
r′
,
z)
dr′
+
a r′
β
(r,
r
′, z − t) dr′];
(3)
n
n
∫
n
n
0
0
∞
1
1
R
(1)
(2)
H r,ϕ,z)
= π(µ-1)
cosnϕ
[
t∫a
(z′)
γ r, 1,
z -tz′
)dz′ + ∫a r′
γ r,
r′,
z)
dr′ +
t
z
∑
n
n
n
n
n=0
0
0
1
1
(3)
(1)
+∫a r′
γ r,
r′,
z-t
)
dr′ +sinnϕ[
t∫b
(z′
)
γ r, 1,
z -tz′)
dz′ +
n
n
n
n
0
0
1
1
(2)
(3)
+∫b r′
γ r,
r′
,
z)
dr′ + ∫b r′
γ r,
r′,
z −t) dr′],
(4)
n
n
n
n
0
0
r
z-d
l
где
r
:
=
,
z
:=
,
t
:=
;
(5)
R
R
R
n+1
2
2
2
(−
1)
bδ a,b,q)
n
(
a
+
b
+
q
)
2ab
2
n
n+1
α
n
(a,b,q):=
[a
−
1
]P ε a,b,q)) +
P
1
(ε(a,b,q));
(6)
1
2
2
2
aΓ(n
−
)
2n
+1
4n
−1
2
n+1
2
2
2
(
−1)
nbδ a,b,q)
a
+
b
+
q
4ab
n
n+1
β
(a,b,q):
=
1
P ε a,b,q))+
P
(ε(a,b,q))
;
(7)
n
1
2
1
aΓ
(n
−
)
2
n
+
1
2
4n
−
1
2
2
n+1
(−
1)
bqδ a,b,q)
n
γ
n
(a,b,q):
=
1
P ε a,b,q)).
(8)
1
2
Γ
(
n
−
)
2
В формулах (6)—(8):
1
2
−
2
−3
2
2
2
2
2
2
4
2ab
δ
(a,b,q)
=(
a
+
b
+
q
)
−
4
a
b
,ε(a,b,q)
=1−
,
(9)
2
2
2
a
+b
+
q
n
где
1
(⋅)
P
— присоединенная функция Лежандра [9], Γ(⋅) — гамма-функция.
2
Дефектоскопия
№ 4
2022
66
В.В. Дякин, О.В. Кудряшова, В.Я. Раевский
∞
(1)
(2)
(3)
Функции {
n
a z),
n
a r),
a
n
(r)
}
в (2) — (4) для каждого n = 0, 1, 2, … есть решение систе-
n=0
мы одномерных линейных интегральных уравнений:
1
1
(1)
(1)
(2)
a z)-2λ π
(z′)α
(1, 1,
t(z- z′))
dz′+
a r′)α
(1,
r′,
tz)
dr′+
n
∫a
n
n
∫
n
n
0
0
1
2
(3)
(01)
+
a r′)α
(1,
r′,
t(1−
z))
dr′
=
a z);
∫
n
n
n
(10)
µ+1
0
1
1
2
(2)
(1)
(3)
(02)
a r)−2λ π
t∫a
(z
′
)γ
(
r, 1,
tz′
)
dz′+
a r′)γ
(r,
r′,
t)
dr′ =
a r);
(11)
n
n
n
∫
n
n
n
µ+1
0
0
1
1
2
(3)
(1)
(2)
(03)
a r)−2λ π
t∫a
(
z
′)γ
(r, 1,
t
(1
−
z′))dz′+
a r′)γ
(r,
r′,
t)dr′ =
a r).
(12)
n
n
n
∫
n
n
n
µ+1
0
0
∞
(1)
(2)
(3)
Функции {
n
b z),
n
b r),
b
n
(r)
}
в (2)—(4) для каждого n = 1, 2, … есть решение системы
n=1
одномерных линейных интегральных уравнений:
1
1
(1)
(1)
(2)
b z)-2λ π
(
z′)α
(1, 1,t(z- z′
))
dz′+
b r′)
α
(1,
r′,
tz)dr′+
n
∫b
n
n
∫
n
n
0
0
1
2
(3)
(01)
+
b r′)α
(1,
r
′
,
t
(1−
z))dr′ =
b z);
(13)
∫
n
n
n
µ+
1
0
1
1
2
(2)
(1)
(3)
(02)
b r)−2λ π
t∫b
(z′
)
γ
(r,1,
tz
′
)dz′+
b r′)γ
(r,
r′,
t)dr′ =
b r);
(14)
n
n
n
∫
n
n
n
µ+1
0
0
1
1
2
(3)
(1)
(2)
(03)
b r)−2λ π
t∫b
(
z
′)γ
(r, 1,
t
(1
−
z′))dz′+
b r′)
γ
(r,
r′ t)
dr′ =
b r).
(15)
n
n
n
∫
n
n
n
µ+1
0
0
В системах (10)—(12) и (13)—(15) параметр λ = (µ - 1)/(µ + 1), а сами системы отличаются
только правыми частями и обозначениями искомых функций. Функции в правых частях этих си-
∞
∞
(01)
(02)
(03)
(01)
(02)
(03)
стем {a z),
a r),
a
(r
)
}
b z),
b r),
b
(r)
суть коэффициенты разложения
n
n
n
и {
n
n
n
}
n
=0
n=1
0
0
0
известных функций
ψ
(ϕ, z),
ψ
(r,ϕ) ,
ψ
(r,ϕ) (связанных с заданным внешним полем H0(r)) в
1
2
3
ряд Фурье вида:
∞
0
(01)
(01)
ψ ϕ,z)
=
a z)cos
nϕ+b z) sinn
ϕ;
(16)
1
∑
n
n
n=0
∞
0
(02)
(02)
ψ ϕ)
=
a r)cos
nϕ+b r) sinn
ϕ;
(17)
2
∑
n
n
n=0
∞
0
(03)
(03)
ψ ϕ)
=
a r)cos
nϕ+b r) sinn
ϕ
,
(18)
3
∑
n
n
n=0
где
0
0
ψ ϕ,z)
=H Rcosϕ,Rsin
ϕ
,
lz + d
),
ϕ∈(
−π π],
z∈
[0, 1];
(19)
1
n
Дефектоскопия
№ 4
2022
Расчет напряженности магнитного поля от однородного цилиндра конечных размеров...
67
0
0
ψ r,ϕ)
=H
(
Rr ϕ
Rrsinϕ,
d),
ϕ∈(
−π π],
r∈[0, 1];
(20)
2
n
0
0
ψ r,ϕ)
=H
(
Rr ϕ
Rrsinϕ,
l+d),
ϕ∈(
−π π],
r∈[0, 1],
(21)
3
n
а
H
0(
x y,
z)
— нормальные составляющие заданного внешнего поля H0(r) = H0(x, y, z) в декарто-
n
вой системе координат на боковой поверхности цилиндра S1 ( в (19)), на его нижнем S2 (в (20)) и
верхнем S3 основании (в (21)), а r = (x, y, z) — декартовы координаты точки.
R
Для расчета декартовых координат поля реакции
H
(r)
в точках
r
=
(0,0, )
z
на оси цилиндра
(внутри или вне его, т.е. z либо внутри интервала (d, d + l), либо вне отрезка [d, d + l]) в работе [4]
получены следующие расчетные формулы:
1
µ-1
R
(1)
2
−32
H
(0,0,z)
=-
(
z
′)[1+ (z -tz′)
]
dz′+
x
∫a
1
4
0
1
1
(2)
2
2
2
−32
(3)
2
2
2
−32
+
a
(r′) r
′
(r
′
+
z
)
dr
a
′ +∫ ′
(r′)⋅ r
[r
′
+
(z -t)
]
dr′;
(22)
∫
1
1
0
0
1
1
µ-1
R
(1)
2
−32
(2)
2
2
2
−32
H
(0,0,z)
=-
(z′)[1+
(z -tz′)
]
dz′+
b r′)r′
(
r′
+
z
)
dr′+
y
∫b
1
∫
1
4
0
0
1
(3)
2
2
2
−32
+
b
(r′)r
′
[
r′
+
(
z -t)
]
dr
′
;
(23)
∫
1
0
1
µ-1
R
(1)
2
−32
H
(0,0,z)
=
(
z
)(
′ ′
z -t⋅z
)[1+(z -tz′)
]
dz′+
z
∫a
0
2
0
1
1
(2)
2
2
−32
(3)
2
2
−32
+
z∫a
(r′)r′(r
′
+
z
)
dr
′ +
(z -t)∫a
(r′)r′
[
r′
+
(z -t)
]
dr′,
(24)
0
0
0
0
(1)
(2)
(3)
(1)
(2)
где функции
{a
(z),a
(r),a
(r)}
— решение системы (10) — (12) для n = 0,
{a
(z),a
(r),
0
0
0
1
1
(3)
(1)
(2)
(3)
a
(r)} — решение системы (10) — (12) для n = 1, а
{b
(z),b
(r),b
(r)} — решение системы
1
1
1
1
(13)—(15) для n = 1.
0
0
0
0
0
0
0
Для случая постоянного внешнего поля
H
=
{H
, H
, H
}
, где
H
, H
, H
— константы,
x
y
z
x
y
z
формулы (2)—(4) для цилиндрических координат поля реакцииHR упрощаются таким образом,
что в (2) и (4) ненулевыми оказываются только два первых члена ряда, а в (3) — только первый
(1)
(2)
(3)
член ряда. При этом участвующие в этих формулах функции
a
(z), a
(r), a
(r) суть решение
0
0
0
системы:
1
1
(1)
(1)
(2)
a
0
(z) - 2λ π
0
(z′)α
0
(1, 1, t(z - z′)) dz′+
a
0
(r′)α
0
(
1, r′,
tz) dr
′+
∫a
∫
0
0
1
(3)
+
a
(r′)α
(
1, r′, t(1− z)
) dr
′
=
0;
(25)
∫
0
0
0
1
1
2
(2)
(1)
(3)
0
a r)−2λ π
t∫a
(z
′
)
γ
(r, 1,
tz′
)
dz
′+
a r′)
γ
(r,
r′,
t
)
dr′
=−
H
;
(26)
0
0
0
∫
0
0
z
µ+
1
0
0
Дефектоскопия
№ 4
2022
68
В.В. Дякин, О.В. Кудряшова, В.Я. Раевский
1
1
2
(3)
(1)
(2)
0
a r)−2λ π
t∫a
(z′)γ
(r, 1,
t(1−
z′))
dz′+
a r′)γ
(r,
r′,
t)
dr′ =
H
,
(27)
0
0
0
∫
0
0
z
µ+1
0
0
(1)
(2)
(3)
где функции
a
(z),a
(r),a
(r)
— решение системы:
1
1
1
1
1
(1)
(1)
(2)
a
(z) - 2λ π
(z′)α
(1, 1, t(z - z′)) dz′+
a
(r′)α
(1, r′, tz) dr′+
1
∫a
1
1
∫
1
1
0
0
1
2
(3)
0
+
a r′)α
(1,
r′,
t(1
−
z))
dr′
=
H
;
∫
1
1
x
(28)
µ+1
0
1
1
(2)
(1)
(3)
a
(r) − 2λ π
t∫a
(z′)
γ
(r, 1, tz′) dz′+
a
(r′)
γ
(r, r′, t) dr′ =
0;
(29)
1
1
1
∫
1
1
0
0
1
1
(3)
(1)
(2)
a
(r) − 2λ π
t∫a
(z′)γ
(r, 1, t(1
− ′)) dz′+
a
(r′)γ
(r, r′, t) dr′ =
0,
(30)
1
1
1
∫
1
1
0
0
(1)
(2)
(3)
а функции
b
(z),b
(r),b
(r) — решение системы:
1
1
1
1
1
(1)
(1)
(2)
b
(z) - 2λ π
(z′)α
(1, 1, t(z - z′)) dz′+
b
(r′)α
(1, r′, tz) dr′+
1
∫b
1
1
∫
1
1
0
0
1
2
(3)
0
+
b r′)α
(1,
r′,
t(1
−
z))
dr′ =
H
;
∫
1
1
y
(31)
µ+1
0
1
1
(2)
(1)
(3)
b
(r) − 2λ π
t∫b
(z′)γ
(r, 1, tz′) dz
′+
b
(r′)γ
(r, r′, t) dr′ =
0;
(32)
1
1
1
∫
1
1
0
0
1
1
(3)
(1)
(2)
b
(r) − 2λ π
t∫b
(z′)γ
(r, 1, t(1
− ′)) dz′+
b
(r′)γ
(r, r′, t) dr′ =
0,
(33)
1
1
1
∫
1
1
0
0
где
bδ b,q)
0
1
α
(1,b,q)
=
P
1
(
ε(1,b,q)
)
−2bP
1
(
ε(1,b,q)
)
;
(34)
0
2
2
2
π
bqδ a,b,q)
0
γ
0
(a,b,q)
=
1
P ε a,b,q));
(35)
2
2
π
2
2
bδ b,q)
2
−
b
−
q
2b
1
2
α
(1,b,q)
=
P
(ε(1,b,q))
+
P
(ε(1,b,q));
(36)
1
1
1
π
3
2
3
2
bqδ a,b,q)
1
γ
1
(a,b,q)
=
1
P ε a,b,q)).
(37)
2
π
Дефектоскопия
№ 4
2022
Расчет напряженности магнитного поля от однородного цилиндра конечных размеров...
69
3. ОСОБЕННОСТИ ПРОГРАММНОЙ РЕАЛИЗАЦИИ
Для реализации указанной в предыдущем параграфе схемы была составлена компьютерная
программа на языке ФОРТРАН. Основным по алгоритмической и аналитической трудоемкости
пунктом этой схемы является нахождение численного решения системы линейных одномерных
интегральных уравнений (10)—(12) (система (13)—(15) отличается от нее только правой частью).
Для дискретизации этой системы интегральных уравнений был применен метод квадратур [10,
с. 157], а в качестве квадратурной формулы, заменяющей интегралы в системе (10)—(12), выбрана
формула средних прямоугольников, которая при разбиении отрезка интегрирования [0, 1] на m
равных частей имеет вид:
1
m
1
f x)dx
≈
f
(x
),
(38)
∫
∑
i
0
m
i=1
где узлы формулы:
2i
−
1
x
=
, i = 1, …, m.
(39)
i
2
m
Прежде чем с помощью (38) и (39) выписать дискретное приближение системы (10)—(12), не-
обходимо преобразовать первый интеграл в левой части (10), поскольку его ядро
(1, 1,
(
))
α
t
z-z′
n
имеет (интегрируемую) логарифмическую особенность при z′ = z. А, именно, можно показать, что
для всех n = 0, 1, … при q→0:
1
α
(1, 1,
q)
=-
(
ln | q | −D
)
+O qln |q|),
n
32
n
4π
(40)
n+1
1
1
D
:= ln8-2
−
-1.
n
∑
k=1
2k
−1
2n
+1
Выделяя стандартным образом указанную особенность [10, с. 159, 160], получим эквивалент-
ный вид уравнения (10):
1
1
(1)
(1)
(1)
(1)
a z)−2λ π
t
a
(z′)
−a z)
α
(1, 1,
t
(z- z′))
dz′+a z)
α
(1, 1,
t(z- z′))
dz′
+
n
∫
(
n
n
)
n
n
∫
n
0
0
1
1
2
(2)
(3)
(01)
+
a r′)α
(1,
r′,
tz)
dr′+
a r′)α
(1,
r′,
t(1−
z))
dr′
=
a z).
(41)
∫
n
n
∫
n
n
n
µ+1
0
0
Заменяя в системе (41), (11), (12) интегралы квадратурными формулами вида (38), (39), после не-
которых преобразований получим для определения дискретных аналогов искомых функций
(1)
(1)
(2)
(2)
(3)
(3)
a
(x
),,
a
(x
);
a
(
x
),,
a
(
x
);
a
(
x
),,a
(x
)
(42)
{
n
1
n
m
n
1
n
m
n
1
n
m
}
следующую систему из 3m уравнений:
m
m
n
(1)
n
(1)
(1)
n
(2)
(q+mp
k
)a
n
(x
k
) +
∑
a
ki
(
a
n
(x
i
)
−
a
n
(x
k
)
)
+
∑
b
ik n
a
(x
i
) +
i=1
i=1
(43)
m
n
(3)
(01)
+
b
a
(
x
)
=
fa
(
x
),
k
=
1,
,
m
;
∑
i,m
+1−k n
i
n
k
i =1
m
m
(2)
n
(1)
n
(3)
(02)
qa
(
x
) +
c
a
(x
)
+
d
a
(
x
)
=
fa
(x
),
k
=
1,
,m;
(44)
n
k
∑
ik n
i
∑
ik n
i
n
k
i
=1
i
=1
Дефектоскопия
№ 4
2022
70
В.В. Дякин, О.В. Кудряшова, В.Я. Раевский
m
m
(3)
n
(1)
n
(2)
(03)
qa
(x
) +
c
a
(x
) +
d
a
(x
)
=
fa
(x
),
k
=1,,m,
(45)
n
k
∑
m+1−i,k n
i
∑
ik n
i
n
k
i
=1
i
=1
которая естественным образом сводится к стандартному виду системы 3m линейных уравнений
относительно 3m неизвестных (42). В формулах (43) — (45) приняты обозначения:
m
0,
k = i
q:=-
,
a
n :
,
k,i =1, ,m;
(46)
ki
=
2λt
π
α
(
1, 1,t(x
−
x
)
)
,
k ≠ i
n
k
i
α
(
1,
x
,
tx
)
γ
(
x
,
x
,
t
)
n
n
i
k
n
n
n
k
i
b
:=
,
c
:=γ
(
x
, 1,
tx
)
,
d
:=
,
k
,
i
=
1,
,
m;
(47)
ik
ik
n
k
i
ik
t
t
1
m
n
f
:=-
,
p
k
:=∫α
n
(
1, 1,
t
(x
k
−x)
)
dx
,
k =1, ,m.
(48)
π
(µ-1)t
0
Интеграл в (48) — сходящийся несобственный 2 рода, поэтому для применения квадратурных
формул он был преобразован с учетом асимптотики (40) к следующему виду:
1
n
1
p
:
=∫ α
(
1, 1,
t
(x
−
x)
)
+
(
lnt+ln |
x
−
x
|
−
D
)
dx
−
k
n
k
32
k
n
4π
0
1
−
[
lnt
−
D
+
x
ln
x
+
(1−
x
) ln(1−
x
)
−
1
]
32
n
k
k
k
k
4π
В этом выражении интеграл особенностей не имеет, поскольку, согласно (40), подынтегральная
функция имеет порядок
O(|
x-x
| ln |
x-x
|)
при
x→x
k
k
k
После аналогичной замены в системе (13) — (15) интегралов квадратурными формулами вида
(38), (39) для определения дискретных аналогов искомых функций
(1)
(1)
(2)
(2)
(3)
(3)
{
b
(x
),
,b
(x
);
b
(
x
),,
b
(
x
);
b
(x
),
,
b
(x
)
}
(49)
n
1
n
m
n
1
n
m
n
1
n
m
получаем ту же систему линейных уравнений (43) — (45) с естественной заменой в правых частях
(01)
(02)
(03)
(01)
(02)
(03)
a
n
(x
k
),
a
n
(x
k
),
a
n
(x
k
)
на
b
n
(x
k
),
b
n
(
x
k
),
b
n
(x
k
)
соответственно.
После решения указанных систем линейных уравнений приближенное вычисление компонент
поля реакции проводится по формулам (2) — (4) и (22) — (24) с естественной заменой интегралов в
них квадратурными формулами вида (38), (39) с использованием полученных дискретных аналогов
искомых функций (42) и (49).
При вычислении матричных элементов системы (43) — (46), а также при последующем вы-
числении компонент поля реакции, используются при различных значениях их параметров функ-
ции
α a,b,q)
,
β a,b,q
)
,
γ a,b,q)
, определенные формулами (6)—(8), каждая из которых со-
n
n
n
n
держит вычисление функции Лежандра
1
(
)
P
z
при z = ε(a,b,q) в (9). Линейные размеры R и l
2
цилиндра могут меняться произвольным образом (для проверки правильности работы програм-
мы рассматриваются и предельные случаи с известными аналитическими ответами бесконечно
длинного цилиндра l→∞ и пластины конечной толщины R→∞), а потому аргумент z функции
Лежандра может принимать сколько угодно большие значения (z→∞) и как угодно близкие к 1
(z→1+0). Для подбора наилучшей вычислительной формулы для этой функции при различных
значениях аргумента проводились численные эксперименты с ее представлением через гиперге-
ометрические функции различного типа (например, [9, с. 128, формула (16)] и [9, с. 130, формула
(24)]), через интегральное представление [9, с. 158, формула (14)], рекуррентные формулы [9,
с. 161, формула (1)]. Наилучшей по точности оказалась комбинированная формула, учитывающая
асимптотическое поведение при z→+∞, а также предельное поведение при z→1 + 0. За основу для
n = 0, 1, 2, … взята формула
n
n
/2
n
(−
1)
z
−
1
1
1
z
−1
1
P z)
=
R
z
+1
Fn
−
,
−
;
n
+
1;
,
(50)
n
2
2
z
+
1
2
2
z
+1
Дефектоскопия
№ 4
2022
Расчет напряженности магнитного поля от однородного цилиндра конечных размеров...
71
n
3
где F(a, b; c; z) — гипергеометрическая функция [9, с. 69]), R0:=1,
R
n
:=
∏
i
−1−
,
n = 1, 2, ….
i=1
4i
При достаточно больших значениях z используется асимптотическая формула
n+1
1
(−1)
2
Γ
n
−
2
2
n+1
2
(4n
−1)lnz
4n
−1
1
1
n
1
P z)
=
z
1−
+
-
3ln 2+2
−
3/2
2
∑
2
π
8z
8
2
2k
−1
k=1
n(n
+1)
1
1
lnz
(51)
−
+
+O
,
2
5/2
2
8
z
z
n
а для значений z, близких к 1, используется соотношение
lim
P
1
(z)
=δ
(δij — символ Кронекера)
n
0
z→ +0
2
[9, с. 164, формула (4)]. Формула (50) получена из представления [9, с. 128, формула (16)]
n/2
n
1
z
+1
1
1
1
z
−1
1
P z)
=
z
+1
F−
,
−
-n
; 1− n;
2
2
z
−1
Γ(1−n)
2
2
z
+1
с учетом соотношения [11, с. 1014, формула (9.101)]
F(αβ;γ;z)
α(α+
1)
(α+ β(
n
β+1)(β+n)
n+1
lim
=
z
F(α+n
+
1,β+
n
+1;n
+2,z).
γ→-n
Γ γ)
(n
+1)!
Формула (51) получается из (50) с использованием соотношения для гипергеометрической
функции, приведенного в [9, с. 117, формула (12)], и некоторых преобразований.
При вычисленииanki в (46) и nkp в (48) возникает необходимость использования формулы (6)
для αn(a, b, q) при a = b = 1 и достаточно малых значениях q. В этом случае особенность типа (40)
вносит достаточную погрешность в численные расчеты по формуле (6). Поэтому для вычисления
αn(1, 1, q) при малых q из (6), (50) и представления гипергеометрической функции в форме, исполь-
зующей ψ-функцию [9, с. 117, формула (12)], в результате достаточно длительных преобразований
получена следующая формула:
n/2
1
1
−
x
B
n
(x)
α
(1,1,q)
=
−
n xqA x)
,
(52)
n
3/2
n
1
(xq)
3/2
2
3/4
4π
(2n+
1)Γ
n
−
(
q
+4)
2
где
2
1
2+
q
x:=
+ ;
2
2
q
4
+
q
∞
1
2n
+
1
2
1
B x):=Γ
n
−
2
n
+1+
1+
c x)
+
1
−
Q x)−Q x);
n
∑
k
n+1
n
2
4
x
k=0
2n
−
1
x
1
2n
-1
1
A x):=8Γ
n
-
1+
-
Q x);
n
2
n
2
4
x
x
∞
2k
+
3
c
(x):=1,
c x):=
c x)
,
k=
0, 1,
;
Q x):
a n,x)
b n,
x);
0
k
+1
k
n
=∑
k
k
2(k
+
3)
x
k=0
Дефектоскопия
№ 4
2022
72
В.В. Дякин, О.В. Кудряшова, В.Я. Раевский
3
3
k
+
k
+
+n
3
2
2
a n,x):=Γ
n
+
,
a n,x):=
a n,x),
k=
0, 1,
;
0
k+1
k
2
(k
+1)(k
+3)x
n+1
1
1
0
b n,x):=-
+ln16
−2
∑
+lnx;
2
i=1
2i
-1
1
1
1
1
b n,x):=
b n,x)+
+
-
-
,
k=
0, 1,
k+1
k
k
+1
k
+3
3
3
k
+
k
+
+n
2
2
Как показали численные эксперименты, формула (52) дает при малых значениях q достаточно
точные результаты для αn(1, 1, q) при замене рядов в Bn(x) и Qn(x) на конечные суммы с относитель-
но небольшим числом слагаемых (в программе хватало 10 слагаемых).
4. ПРОВЕРКА НА ИЗВЕСТНЫХ ПРЕДЕЛЬНЫХ СЛУЧАЯХ
В виду сложности и объемности аналитических преобразований, приведших к описанному выше
алгоритму вычисления напряженности поля реакции (а потому и результирующего поля) от однород-
ного цилиндра, а также объемности реализующей этот алгоритм компьютерной программы, важное
значение приобретает качественная и количественная проверка результатов работы этой программы.
Было проделано качественное тестирование, предполагающее проверку наличия определенной
симметрии вычисляемой напряженности результирующего поля в зависимости от типа симметрии
напряженности поля внешнего источника. В частности, если напряженность постоянного внешне-
го поля H0 направлена по оси z, то для компонент напряженности результирующего поля в любой
точке вне или внутри цилиндра должно быть выполнено Hφ = 0, а Hz и Hr не зависят от φ. Если
напряженность такого поля H0 направлена по оси х, то при смене координаты φ точки наблюдения
на (-φ) декартовы координаты Hx и Hz не меняются, а Hy меняет знак. Если же напряженность H0
направлена по оси у, то при смене координаты φ точки наблюдения на π-φ координаты Hy и Hz не
меняются, а Hx меняет знак. Кроме того, при значительном удалении точки наблюдения от магне-
тика результирующее поле должно практически совпадать с полем внешнего источника.
Было проведено и многократное количественное тестирование результатов работы программы,
заключающееся в количественном сравнении их в предельных случаях цилиндрической формы
магнетика (бесконечно длинный цилиндр, бесконечная пластина конечной толщины), для которых
в случае погружения в постоянное внешнее поле известны точные аналитические формулы напря-
женности результирующего поля вне и внутри магнетика. Ниже приводятся типичные результаты
0
0
0
0
такого сравнения. Напряженность внешнего поля
H
={H
, H
, H
}
при указанном тестировании
x
y
z
бралась постоянной как внутри магнетика, так и в той его окрестности, где производится вычис-
0
0
0
ление напряженности результирующего поля:
H
=10,
H
=
20,
H
=
30.
Цилиндр рассматри-
x
y
z
вался симметричным относительно плоскости z = 0 (то есть параметр d = -l / 2).
Для бесконечно длинного цилиндра с осью вдоль оси z использовались следующие формулы
для напряженности H = {Hx, Hy, Hz} постоянного результирующего поля внутри магнетика [12, с.
209], [13, с. 339]:
0
0
2H
2
H
x
y
0
H
=
,
H
=
,
H
=H
(53)
x
y
z
z
1+µ
1
+µ
Программе задавались следующие параметры: μ = 19, R = 2, l = 250, m = 800, цилиндрические
координаты точки наблюдения r = 1,5, φ = π/3, z = 3. По формуле (53) получаются значения Hx = 1,
Hy = 2, Hz = 30; программа дала Hx = 0,9994, Hy = 1,999, Hz = 29,93.
Напряженность результирующего поля вне указанного бесконечно длинного цилиндра в точке
наблюдения с декартовыми координатами (x, y, z) [14, с. 245]:
0
0
2
2
0
0
0
2
2
0
H
=H
+pH
(x
−y
) + 2H
xy,
H
=H
+pH
(y
−x
) + 2
H
xy,
H
=H
0,
(54)
x
x
x
y
y
y
y
x
z
z
Дефектоскопия
№ 4
2022
Расчет напряженности магнитного поля от однородного цилиндра конечных размеров...
73
2
2
2
2
где
p := λR
/(x
+y
)
. Программе задавались те же параметры, а точка наблюдения вне цилин-
дра имеет координаты r = 3, φ = π/3, z = 3. Формула (54) дает результаты Hx = 14,928, Hy = 27,464,
Hz = 30; программа — Hx = 14,930, Hy = 27,469, Hz = 29,925.
Рассмотрим однородный магнетик, занимающий область Ω между двумя параллельными пло-
скостями с уравнениями z = -l/2 и z = l/2 (бесконечная пластина конечной толщины) . Форма такой
пластины является предельным случаем (радиус основания R →∞) соответствующего цилиндра
высоты l. Напряженность результирующего поля внутри и вне пластины [13, с. 339]:
1
0
0
0
0
0
H
=H
0,
,
H
=H
0,
H
= H
,
r
∈Ω;
H=H
={H
, H
, H
},
r∉Ω
(55)
x
x
y
y
z
z
x
y
z
µ
Программе были заданы следующие параметры: μ = 20, R = 800, l = 6, m = 800; цилиндрические
координаты точки наблюдения внутри цилиндра r = 3, φ = π/3, z = 2, а вне цилиндра r = 3, φ = π/3,
z = 5. По формулам (55) значения напряженности результирующего поля внутри цилиндра Hx = 10,
Hy = 20, Hz = 1,5, а вне его Hx = 10, Hy = 20, Hz = 30. Программа дала соответственно Hx = 9,65,
Hy = 19,3, Hz = 1,54 и Hx = 9,65, Hy = 19,3, Hz = 30,12.
Следующий проделанный вид тестирования работы программы, реализующей предложен-
ный алгоритм, касается проверки выполнения необходимых граничных условий. В точках гра-
ницы S магнетика для нормальной составляющей Hn напряженности результирующего поля
(i)
(e)
должно выполнятьсяH(e)
= µH
, где
H
и H(i)n суть предельные значения Hn на поверхности
n
n
n
S извне и изнутри магнетика соответственно [13]. Таким образом, должно выполняться:
( )e
H
n
=µ
(56)
(i
)
Hn
Проверка программы в отношении выполнения (56) проводилась для симметричного относи-
тельно плоскости z = 0 цилиндра с параметрами μ = 19, радиус основания R = 2, высота l = 6,
узлов квадратурной формулы m = 800.
Для боковой поверхности цилиндра S1 нормальная составляющая совпадает с r-компонентой
напряженности, поэтому (56) имеет вид:
( )
e
H
r
=µ
(57)
(i
)
Hr
Для проверки (57) составляющаяH(e)r вычислялась в точке с цилиндрическими координатами
r = 2,01, φ = π/3, z = 2, аH(i)r — в точке r = 1,99, φ = π/3, z = 2. Полученное программой отношение
компонент в (57) равно 18,86, то есть погрешность составила 0,7 %.
Для верхнего основания цилиндра S3 нормальная составляющая совпадает с z-компонентой
напряженности, поэтому (56) имеет вид
( )
e
H
z
=µ
(58)
(i
)
Hz
Для проверки (58) составляющаяH(e)z вычислялась в точке с цилиндрическими координатами
r = 1,5, φ = π/3, z = 3,01, аH(i)z — в точке r = 1,5, φ = π/3, z = 2,99. Полученное программой отно-
шение компонент в (58) равно 18,66, то есть погрешность составила 1,8 %.
Проведенные проверки подтверждают справедливость предлагаемого алгоритма вычисления
результирующей напряженности поля и его программной реализации. Отметим, что при сравне-
нии эффективности работы описанной компьютерной программы с программой, реализующей
алгоритм в работе [3], оказалось, что первая дает значительно более точные результаты за значи-
тельно более короткое время.
5. ЗАКЛЮЧЕНИЕ
Кратко сформулируем основные результаты, полученные в настоящей работе.
1. Приведен алгоритм нахождения напряженности результирующего магнитного поля внутри
и вне однородного цилиндра конечных размеров, помещенного во внешнее магнитное поле про-
извольной конфигурации, сводящийся в основной его части к решению некоторого количества си-
стем трех одномерных линейных интегральных уравнений.
Дефектоскопия
№ 4
2022
74
В.В. Дякин, О.В. Кудряшова, В.Я. Раевский
2. Описаны подводные камни и способы решения возникающих проблем при составлении и
программной реализации приведенного алгоритма дискретизации указанных систем уравнений.
В частности, уравнения системы приведены к виду, в котором отсутствует изначальная сингуляр-
ность в ядрах интегральных операторов, выведены подходящие асимптотические формулы для вы-
числения функций Лежандра и выражений, содержащих такие функции и входящих в ядра инте-
гральных операторов, что позволило повысить точность получаемых результатов и применимость
описанного алгоритма для всего многообразия возможных геометрических параметров цилиндри-
ческого магнетика.
3. Приведены данные качественного и количественного тестирования результатов работы реа-
лизующей описанный алгоритм компьютерной программы, состоящего в проверке их соответствия
физическим законам явления, а также в количественном сравнении этих результатов в предельных
частных случаях с известными аналитическими ответами.
4. Проведено сравнение эффективности работы компьютерных программ, реализующих под-
ход настоящей статьи и подход, описанный в нашей предыдущей работе [3]. Основные выводы
такого сравнения следующие:
расчет с приемлемой точностью напряженности результирующего поля в рамках подхода на-
стоящей статьи требует на стандартном бытовом персональном компьютере от нескольких секунд
до 40 секунд (в зависимости от задаваемого числа разбиений в используемых квадратурных форму-
лах), в то время как аналогичный расчет в рамках подхода [3] занимает от 5 до 10 мин;
расчет в рамках подхода настоящей статьи позволяет получать напряженность результирую-
щего поля с хорошей точностью даже в точках непосредственной близости к границе магнетика, а
расчет в рамках подхода [3] при приближении точки наблюдения к границе магнетика резко теряет
в точности.
Работа выполнена в рамках государственного задания по теме
«Квант»
(“Quantum”)
№ АААА-А18-118020190095-4.
СПИСОК ЛИТЕРАТУРЫ
1. Сапожников А.Б. Теоретические основы магнитной дефектоскопии металлических тел. Томск:
Изд-во ТГУ, 1980. 308 с.
2. Дякин В.В., Кудряшова О.В. Дефект в цилиндре // Дефектоскопия. 2012. № 4. С. 41—55.
3. Дякин В.В., Кудряшова О.В., Раевский В.Я. К расчету поля конечного магнитного цилиндра //
Дефектоскопия. 2019. № 10. С. 24—34.
4. Dyakin V.V., Kudryashova O.V. , Raevskii V.Ya. One Approach to the Numerical Solution of the
Basic Equation of Magnetostatics for a Finite Cylinder in an Arbitrary External Field // Russian Journal of
Nondestructive Testing. 2021. V. 57. No. 4. P. 291—302. [Дякин В.В., Кудряшова О.В., Раевский В.Я. Один
подход к численному решению основного уравнения магнитостатики для конечного цилиндра в произ-
вольном внешнем поле // Дефектоскопия. 2021. № 4. С. 22 — 34.]
5. Дякин В.В., Кудряшова О.В., Раевский В.Я. О проблемах использования пакетов универсальных
программ для решения задач магнитостатики // Дефектоскопия. 2018. № 11. С. 23—34.
6. Хижняк Н.А. Интегральные уравнения макроскопической электродинамики. Киев: Наукова дум-
ка, 1986. 280 с.
7. Дякин В.В. Математические основы классической магнитостатики. Екатеринбург: РИО УрО РАН,
2016. 404 с.
8. Friedman M.J. Mathematical study of the nonlinear singular integral magnetic field equation. 1. // SIAM
Journal Appl. Math. 1980. V. 39. № 1. P. 14—20.
9. Бейтмен Г., Эрдейи А. Высшие трансцендентные функции. Т. 1. М.: Наука, 1973. 294 с.
10. Верлань А.Ф., Сизиков В.С. Интегральные уравнения: Методы, алгоритмы, программы / Спра-
вочное пособие. Киев: Наукова думка, 1986. 544 с.
11. Градштейн И.С., Рыжик И.М. Таблицы интегралов, рядов, произведений. Санкт-Петербург:
БВХ-Петербург, 2011. 1232 с.
12. Татур Т.А. Основы теории электромагнитного поля. М.: Высшая школа, 1989. 271 с.
13. Ахиезер А.И. Общая физика. Электрические и магнитные явления. Киев: Наукова думка, 1981.
471 с.
14. Неразрушающий контроль и диагностика / Под ред. В.В. Клюева М.: Машиностроение, 1995.
487 с.
Дефектоскопия
№ 4
2022
