Электромагнитные методы
УДК 620.179.14
К РАСЧЕТУ ПОЛЯ КОНЕЧНОГО МАГНИТНОГО ЦИЛИНДРА
© 2019 г. В.В. Дякин1, О.В. Кудряшова1,*, В.Я Раевский1,**
1 Институт физики металлов имени М.Н. Михеева Уральского отделения Российской академии наук,
Россия 620108 Екатеринбург, ул. С. Ковалевской, 18
E-mail: *kudryashova_ov@imp.uran.ru**raevskii@imp.uran.ru
Поступила в редакцию 08.04.2019; после доработки 17.05.2019
Принята к публикации 24.05.2019
Разработан алгоритм численного решения прямой линейной задачи магнитостатики по вычислению поля от конеч-
ного цилиндра с постоянной магнитной проницаемостью μ. В ходе численных экспериментов установлены параметры
дискретизации системы уравнений, выполнены проверки решения на соответствие результатам точно решаемых задач,
обозначены границы применимости бесконечной модели цилиндра и намечены перспективы дальнейшего использова-
ния численного алгоритма и компьютерной программы.
Ключевые слова: теоретические задачи магнитостатики, магнитная дефектоскопия, интегродифференциальное уравне-
ние магнитостатики, линейная задача магнитостатики, поле магнитного цилиндра, реализация метода конечных элементов.
DOI: 10.1134/S0130308219100038
ВВЕДЕНИЕ
В научном обороте широко используются бесконечные модели, призванные служить неким
приближением при рассмотрении в чем-то похожих конечных тел. При этом делается ряд огово-
рок, обосновывающих такую замену и убеждающих считать ее оправданной. И действительно,
бесконечные модели приносят немалую пользу, например, могут упростить задачу, сделать ее раз-
решимой аналитически или хотя бы привести к такому виду, на основе которого может быть
построена эффективная по точности решения и приемлемым затратам машинного времени и
машинной памяти вычислительная процедура.
Что касается тел с цилиндрической симметрией, хорошо известна задача о бесконечном магнит-
ном цилиндре с постоянной магнитной проницаемостью, помещенном в бесконечную магнитную
среду с иной магнитной проницаемостью в однородном вешнем поле, перпендикулярном оси цилин-
дра. Решение здесь представляется в элементарных функциях [1, с. 245]. Это самая первая теорети-
ческая задача магнитостатики, которая была приспособлена для нужд магнитной дефектоскопии еще
в 30-х годах прошлого столетия. В [2] можно найти решение по расчету поля бесконечного магнит-
ного цилиндра при условии неоднородного намагничивания. Выражения для компонент результиру-
ющего поля уже не является столь простыми и записываются через специальные функции.
Осуществлялись попытки и по внедрению дефекта вовнутрь цилиндра. Ранее мы обращались
к решению задачи линейной магнитостатики для тела в виде цилиндра, содержащего внутреннее
несоосное включение такой же цилиндрической формы, но отличающееся значением магнитной
проницаемости. Используя переход к бицилиндрическим координатам, было получено аналитиче-
ское решение [3, 4]. Кроме того, к решению данной задачи был применен и один полуаналитиче-
ский метод, позволяющий в реальном времени получать результат с контролируемой точностью
[5]. Второй подход был использован и для решения задачи о полом цилиндре (трубе) с бесконечно
протяженным цилиндрическим дефектом [6]. В [3—6] дефектная область также является неогра-
ниченной, как и то тело, в котором она находится.
Естественно возникает потребность рассмотреть задачу магнитостатики для следующих тел:
однородный конечный цилиндр (а), конечный цилиндр с внутренним (разумеется, конечным) дефек-
том (б), однородный полый конечный цилиндр (труба) (в), конечный фрагмент трубы с внутренним
дефектом (г). Выполнение задач именно в этом порядке позволит тестировать последующую на
основе предыдущей и даст возможность оценить качество и точность результата — предстоит осу-
ществить численное решение. Сейчас мы приступаем к работе над первым пунктом программы.
Для реализации этой задачи мы обратились к методу конечных элементов [7]. Предположение
о постоянстве магнитной проницаемости позволило свести трехмерную задачу к двумерной и тем
самым обеспечить практически приемлемые точность решения и время расчета.
Рассматриваемая в данной работе задача не только может быть использована как инструмент
тестирования более сложных перечисленных выше задач, но и как инструмент для решения обрат-
К расчету поля конечного магнитного цилиндра
25
z
ной задачи по определению на основе анализа зависимостей
компонент магнитного поля от координат геометрических и
n3
физических параметров магнитного цилиндра, к которому нет
h2
S3
доступа.
1. ПОСТАНОВКА ЗАДАЧИ
Ω: μ > 1
Пусть имеется тело Ω, ограниченное поверхностью S, гео-
0
y
метрический образ которой представлен на рис.1.
Цилиндрический фрагмент поверхности обозначим через S1
ρ0
(радиус цилиндра равен ρ0), а круги в плоскостях z = -h1 и z = h2
n
x
1
S1
обозначим через S2 и S3 соответственно. Примем, что тело Ω
является магнетиком с магнитной проницаемостью μ, не зави-
сящей ни от величины поля, ни от координат: μ = const.
–h1
S
2
Считаем, что магнетик помещен во внешнее намагничивающее
n
поле H0(r), вид которого пока не оговариваем. Поставим задачу
2
R3\Ω: μ = 1
найти результирующее поле.
В качестве математической основы возьмем интегродиффе-
Рис. 1.
ренциальное уравнение магнитостатики [8], записанное относи-
тельно вектор-функции напряженности магнитного поля H(r),
которое для магнетика с постоянной магнитной проницаемостью μ, ограниченного поверхностью
S, имеет вид
µ-1
Hr')
n
0
3
H r)+
dS'
=
H
(
r
)
,
r∈
S,
(1)
4π
r-r'
S
где
(
)
H
r'
составляющая вектора H, взятая на поверхности S изнутри области Ω (на это ука-
n
зывает верхний символ в обозначении) вдоль внешнего по отношению к этой области единичного
вектора нормали n, построенного в точке r S.
Вид уравнения (1) определяет следующую последовательность действий: сначала нужно найти
значение функцииHn на S, после чего выполнить подстановкуHn в (1). При этом (1) становит-
ся выражением по вычислению поля в произвольной точке пространства за исключением точек,
принадлежащих S.
Для вывода уравнения относительно функции
(
)
H
r
, r S, умножим обе части (1) на вектор
n
n и положим, что r Ω. Получим уравнение
µ-1
H
(r)
n
0
H
n
(r)
+
dS′= H
n
(
r
)
,
(2)
4π ∂n
r
r
в котором осуществим предельный переход при r S, используя формулу предельного значения
для нормальной производной потенциала простого слоя [9]
rS
H
n
(r)
1
dS
′ →
2πH
(r)
+
H
(
r
)
dS
(3)
n
n
n
r r
n
r r
S
S
Таким образом, из (2) получаем окончательное уравнение относительно значения нормальной
составляющей искомого поля на внутренней стороне поверхности S
χ
1
2
0
H
n
(r)
+
H
n
(r)
dS′=
H
n
(
r
)
,
rS,
(4)
2
π
n
r
r
µ+
1
S
где χ := (μ - 1)/(μ + 1).
Учитывая, что поверхность S есть объединение поверхностей S1, S2 и S3, на основе (4) записы-
ваем систему интегральных уравнений относительно нормальных составляющих на каждом из
этих фрагментов
Дефектоскопия
№ 10
2019
26
В.В. Дякин, О.В. Кудряшова, В.Я. Раевский
3
χ
1
2
0
H
(r)
+
H
(r)
dS′ =
H
(
r
)
,
rS
,
i
=1, 2, 3,
(5)
n
i
∑∫
n
j
j
n
i
i
2π
n
r
r
µ+1
j=1 S
j
i
где ni (i = 1, 2, 3) — единичный вектор нормали на поверхности с соответствующим номером.
2. ДИСКРЕТИЗАЦИЯ СИСТЕМЫ И ЕЕ ЧИСЛЕННОЕ РЕШЕНИЕ
Решение системы (5) будем искать численными методами. Для этого построим разбиения:
N
h
+h
z
1
2
{z }
:
z
=-h
;
z
=
z
+∆
;
=
— шаг; Nz — число точек, Nz ≥ 2;
l
l=1
1
1
l+1
l
z
z
N
1
z
Nϕ
2π
{
ϕ
}
:
ϕ
=
0;
ϕ
+∆
;
=
— шаг; Nφ — число точек, Nφ ≥ 2;
l
l=1
1
l+1
l
ϕ
ϕ
N
-1
ϕ
Nρ
ρ
0
{
ρ
}
:
ρ
=
0;
ρ
+∆
;
=
— шаг; Nρ — число точек, Nρ ≥ 2.
l
l=1
1
l+1
l
ρ
ρ
Nρ
1
Таким образом, цилиндрическая поверхность оказывается разбитой на совокупность элементов
{
S
}N1
,
где N1 = (Nz - 1)(Nφ - 1), а нижняя и верхняя грани — на множества элементов {S
}N2
,
и
1m m=1
2m m=1
{
S
3
}N3
,
где N2 = N3 = (Nρ - 1)(Nφ - 1); N := N1 + 2N2 — общее число фрагментов. Через rim
m m=1
обозначим радиус-вектор «центра» отдельно взятой m-й площадки на i-м элементе поверхности,
цилиндрические координаты которого принимаем как среднее арифметическое соответствую-
щих граничных значений координат на этой площадке. Далее предположим, что значение нор-
мальной составляющей вектора напряженности магнитного поля является неизменным внутри
отдельно взятой площадки:
H
n
(r)
H
n
(
r
im
), если r Sim, m = 1, 2, …, Ni, i = 1, 2, 3.
i
i
В результате дискретизации система интегральных уравнений (5) сводится к системе линей-
ных алгебраических уравнений относительно неизвестных
H
n
(
r
im
)
i
3
Nj
χ
2
(
i,
j
)
0
H
r
+
H
r
=
H
r
,
(6)
n
(
im
)
∑∑
mk
n
(
jk
)
n
(
im
)
i
j
i
2π
j
=1
k
=1
µ
+
1
где(i,j)mk
— матричные элементы системы:
(i,
j
)
1
:
=
dS
,
(7)
mk
jk
n
r r
S
i
jk
r=r
im
и где на
S
:=
=
;
на
S
:=
=-
;
на
S
:=
=
1
2
3
n
∂ρ
n
z
n
z
1
2
3
Матрица коэффициентов системы состоит из девяти блоков
N
;
N
N
;
N
N
;
N
(
1,1
)
1
1
(
1,2
)
1
2
(
1,3
)
1
3
{
mk
}
{
mk
}
{
mk
}
m=1;k
=1
m=1;k=1
m
=1;k=1
N
2
;
N
1
N
2
;
N
2
N
2
;
N
3
(
2,1
)
(
2,2
)
(
2,3
)
(8)
{
mk
}
{
mk
}
{
mk
}
m
=
1;k
=1
m
=
1;k=1
m
=
1;k=1
N
3
;
N
1
N
3
;
N
2
3
N ; N
3
(
3,1
)
(
3,2
)
(
3,3
)
{
mk
}
{
mk
}
{
mk
}
m=
1;k
=1
m=
1;k
=1
m=
1;k=1
N
1
;
N
1
(
1,1
)
1). Запишем выражения матричных элементов блока
{
mk
}
,
при вычислении которых
m=1;k=1
имеем в виду, что как точки r1m
= (ρ0, φ1m, z1m), так и область интегрирования
Дефектоскопия
№ 10
2019
К расчету поля конечного магнитного цилиндра
27
ϕ
z
S
=
r′=
(
ρ
,ϕ′, z
)
:ϕ′∈ϕ
−δ
,ϕ
и
z
[
z
−δ
, z
]
,
δ
=
1k
{
0
1k
ϕ
1k
ϕ
1k
z
1k
z
}
где
ϕ
и
δ
z
=
, нахо-
2
2
дятся на цилиндрической поверхности. При m k имеем
ϕ
1k
ϕ
z
1k
z
ϕ
−ϕ′
dz
(1,1
)
2
2
1m

=-2ρ
dϕ′sin
=
mk
0
3
2
ϕ
1k
−δ
ϕ
z
1k
−δ
z
ϕ
−ϕ′
2
2
2
2
1m
4ρ
sin
+
(
z
z
)
0
1m
2
t=z
1k
z
1m z
t
+
=-
F
ψ
,
κ
F
ψ
,κ
,
(
mk
(t))
(
mk
(t))
2
2
t=z
z
−δ
4
ρ
+t
1k
1m z
0
2
ψ
4ρ
π-ϕ
±
δ
d
ϑ
2
0
±
1m
1k
ϕ
где
κ
(t)
:=
;
ψ
:=
;
F
(
ψ
,
κ
)
=
— неполный эллиптиче-
2
2
mk
2
2
4ρ
+t
2
0
0
1
sin
ϑ
ский интеграл первого рода.
При m = k интеграл является несобственным, но сходящимся. Для его вычисления окружаем
точку r1m окрестностью, совпадающей по форме с элементарной площадкой, вычисляем возника-
ющие интегралы и переходим к пределу, стягивая окрестность в точку r1m. Таким образом, полу-
чаем выражение
2δ
π+δ
π-δ

(1,1
)
z
ϕ
ϕ
=-
F
,
κ-
F
,
κ
,
mm
2
2

(9)
4ρ
2
2
0
z
2
4ρ
2
0
где
κ
:=
2
2
4ρ
0
z
N
i
;
N
1
(
i,1
)
2). Для матричных элементов блоков
{
mk
}
,
где i = 2, 3, точка rim= (ρim, φim, ηi) принад-
m=1;
k=1
лежит нижней (где i = 2 и η2 = -h1) или верхней грани (где i = 3 и η3 = h2), а интегрирование
осуществляется по фрагменту цилиндрической поверхности, поэтому окончательное выражение
также легко записывается через неполные эллиптические интегралы первого рода:
i
2ρ
(
1
)
z′=z
1k
z
(
i,1
)
0
+
=
F
(
ψ
,κ
(z))
F
(
ψ
,
κ
(
z
))
,
mk
mk
mk
2
2
2
z′=z
−δ
1k
z
ρ
+
(
z′-η
)
+
2
ρ
ρ
0
im
i
0
im
4ρ
ρ
π-ϕ
±
δ
2
0
im
±
im
1k
ϕ
где
κ
(
z
)
:
=
и
ψ
:
=
2
2
mk
ρ
+
(
z′-η
)2
+2
ρ
ρ
2
0
im
i
0
im
N
1
;
N
j
(
1,
j
)
3). При записи расчетных формул матричных элементах блоков
{
mk
}
,
где j = 2, 3, учи-
m=1;k=1
тываем, что точка r1m = (ρ0, φ1m, z1m)S1mS1, а интегрирование осуществляется по фрагменту Sjk
ρ
плоской грани:
S
=
r
′ =
ρ′
,
ϕ′
,η
:
ρ′∈ ρ
−δ
,
ρ
иϕ′∈ϕ
−δ
,
ϕ
,
где
δ
=
,
jk
{
(
j
)
jk
ρ
jk
ρ
jk
ϕ
jk
ϕ
}
ρ
2
η2 = -h1 и η3 = h2. Здесь легко снимается интегрирование по переменной ρ [10, с. 97]:
ϕ
jk
ϕ
ρ
jk
ρ
2
ρ
ρ′
−ρ′
cos
(
ϕ
−ϕ
)
(
1,
j
)
0
1m

=
d
ϕ′
dρ′
=
mk
3
ϕ
jk
−δ
ϕ
ρ
jk
−δ
ρ
2
2
2
2
ρ
+ρ′
+
(
z
−η
)
2
ρ
ρ′
cos
(
ϕ
− ϕ′
)
(
0
1
m
j
0
1m
)
ϕ
jk
ϕ
ρ′=ρ
jk
ρ
=
dϕ′F
(
ρ′,
ϕ′
;
ρ
,
ϕ
,
z
,
η
)
,
0
1
m
1m j
ρ′=ρ
jk
−δ
ρ
ϕ
jk
−δ
ϕ
Дефектоскопия
№ 10
2019
28
В.В. Дякин, О.В. Кудряшова, В.Я. Раевский
где
F
(
ρ′,
ϕ′ρ,ϕ,z,η
)
:=
2
2
2
2
ρsin
ϕ-ϕ′
ρ
+
z
−η
2ρρ′cos
ϕ-ϕ′
−ρ′
z
−η
cos
ϕ-ϕ′
(
)
(
(
)
(
)
)
(
)
(
)
=
+
2
2
2
2
2
2
(
ρ
sin
(
ϕ-ϕ′
)
+
(
z
−η
)
)
ρ
+
ρ′
+
(
z
−η
)
2ρρ′cos
(
ϕ-ϕ′
)
2
2
2
+cos
(
ϕ-ϕ′
)
ln
2
ρ
+ρ′
+
(
z
−η
)
2ρρ′
cos
(
ϕ-ϕ′
)
+2ρ′-2ρcos
(
ϕ-
ϕ′
)
(10)
Подынтегральная функция уже не так проста, однако полученный определенный интеграл
тоже может быть представлен в виде комбинации неполных эллиптических интегралов, но уже и
первого, и второго, и третьего родов. Выражение выглядит крайне громоздким, мы не будем при-
водить его здесь.
N ; N
N
;
N
(
2,2
)
2
2
(
3,3
)
3
3
4). Блоки
и
полностью состоят из нулевых элементов.
{
mk
}
{
mk
}
m=
1;
k=1
m=1;k=1
N
j
(
i,
j
)
5). Матричные элементы блоков
{
mk
}Ni;
, где (i, j) = (2, 3) или (3, 2), определяются из
m=
1;k=1
выражений:
ϕ
jk
ϕ
ρ
jk
ρ
(
i,
j
)
ρ′dρ′
=-
(
h
+h
)
d
ϕ
=
mk
1
2
3
ϕ
ρ
2
2
2
jk
ϕ
jk
ρ
2
ρ
+ ρ′
+
(
h
+h
)
2ρ
ρ′cos
(
ϕ
−ϕ
)
(
im
1
2
im
im
)
ϕ
jk
ϕ
ρ
′=ρ
jk
ρ
=
(
h
+
h
)
dϕ′G
ρ
,ϕ′;ρ
,ϕ
,η
,η
=
1
2
(
im
im i
j
)
ρ′
jk
−δ
ρ
ϕ
jk
−δ
ϕ
ρ′=ρ
jk
ρ
h +h
1
ρ
′ 
ψ
+
mk
1
ρ
ψ
+
mk
1
2
=-
ν
1
+
-
Π
(
ψ
,
ν
,
κ
)
−ν
1+
+
(
ψ
,ν
,
κ
)
,
1;im
1;im
im
2;
im
Π
2;im
im
-
ν
ρ
ψ
mk
ν
ρ
ψ
mk
2
a
+b
im
im
im
im
im
im
ρ′=ρ
−δ
jk
ρ
2
2
ρ
+
(
z
−η
)
− ρρ′cos
(
ϕ-ϕ′
)
где
G
(
ρ′,
ϕ ρ,ϕ,z,
η
)
:
=
;
(11)
2
2
2
2
2
2
(
ρ
sin
(
ϕ-ϕ′
)
+
(
z
−η
)
)
ρ
+ρ′
+
(
z
−η
)
2ρρ′
cos
(
ϕ-ϕ′
)
2
2
2
2
b
im
a
im
:
im
+ρ′
+
(
h
1
+
h
2
)2
;
b
im
:
=
2ρ
im
ρ′;
κ
im
:
=
;
a
+b
im
im
2
ρ
2
2
im
ν
:
=
;
ν
:
=-
;
ν
:
=
;
im
1;im
2;im
(
h
+
h
)2
1
1
1
2
1
+
1
+
1
1+
νim
ν
im
π-ϕ
±δ
±
im
jk
ϕ
ψ
:
=
;
mk
2
ψ
d
ψ
Π
(
ψ ν
κ
)
:
=
— неполный эллиптический интеграл третьего рода.
2
2
2
0
(
1
sin
ψ
)
1−κ
sin
ψ
Элементы вектора-столбца правой части системы (6) вычисляются на основе заданного выра-
0
0
жения вектора напряженности внешнего намагничивающего поля H0(r):
H
(
r
)
=
(
H
(
r
)
,
n
)
,
n
i
im
im
im
где nim — единичный вектор внешней к области Ω нормали, построенной в
точке с радиус-векто-
ром rim, принадлежащей элементу разбиения Sim, i = 1, 2, 3, m = 1, 2, …, Ni. Если подробнее, то n1m=
= (cos φ1m, sin φ1m, 0), n2m = (0, 0, -1) и n3m = (0, 0, 1). При расчетах мы остановились на случаях
однородного внешнего поля с абсолютной величиной H0, ориентированного либо вдоль оси ox,
либо вдоль оси oz. Численные значения компонент результирующего поля выдаются нормирован-
0
ными H0 и вводится обозначение
H=H
/H
На языке программирования Fortran с использованием стандартной программы решения систе-
мы линейных алгебраических уравнений с оценкой числа обусловленности C [11] были составле-
Дефектоскопия
№ 10
2019
К расчету поля конечного магнитного цилиндра
29
ны программы решения системы (6) со следующими формами представления матричных элемен-
N
j
(
i,
j)
тов
(i = 1, 2, 3, и j = 1, 2, 3): (1) — из расчета по формуле приближенного значения
{
}Ni;
mk
m=1;k=1
двойного интеграла как произведения значения подынтегральной функции в точке rjk Sjk на
площадь элемента Sjk (j = 1, 2, 3 и k = 1, …, Nj); (2) — из расчета по формуле приближения соот-
ветствующего определенного интеграла; (3) — из расчета через неполные эллиптические интегра-
N
1
(
1,1
)
лы. При этом во всех случаях матричные элементы вида
,
расположенные на диагонали
{
}
mm
m=1
данного блока, вычислялись исключительно по формуле (9). В ходе расчетов использовались
параметры разбиения, приведенные в табл. 1.
Таблица
1
Вариант
Nz
Nφ
Nρ
N1
N2
N
C <
Время счета <
1
4
6
3
24
18
60
100
1 с
2
10
10
5
100
50
200
100
2 с
3
20
20
10
400
200
800
200
20 с
4
25
25
15
625
375
1375
300
1 мин
5
30
30
20
900
600
2100
300
2 мин
6
40
40
20
1600
800
3200
300
5 мин
7
50
50
20
2500
1000
4500
300
12 мин
8
60
60
20
3600
1200
6000
300
25 мин
9
70
70
30
4900
2100
9100
400
80 мин
При выполнении численных экспериментов установлено, что для всех форм представления
матричных элементов, начиная с шестого варианта происходит стабилизация двух значащих цифр
окончательного значения результирующего поля (формулы для расчета его компонент приводятся
в следующем пункте) при внешнем поле, направленном вдоль оси цилиндра, и четырех — при
внешнем поле, направленном поперек этой оси. Выполненные численные расчеты отразили оче-
видное: с ростом числа элементов разбиения приближение как двойного, так и определенного
интеграла по формуле среднего значения становится более точным. Таким образом, оказывается
достаточным использовать первый или второй варианты представления матричных элементов. Мы
выбрали второй из них, для этого же случая в последней колонке таблицы указано время расчета
на ПК Pentium с достаточно скромными возможностями: тактовой частотой 2,67 ГГц и оператив-
ной памятью 1 ГБ. Было бы весьма интересно ознакомится с результатами решения аналогичной
задачи, полученными с помощью пакета стандартных программ Ansys, а также сравнить затраты
машинных ресурсов. К сожалению, не удалось найти публикации на эту тему.
результирующего поля на площадках разбиения.
На рис. 2 представлены графики зависимости
H
n
(
z
)
для всех вариантов расчета, указанных в
табл. 1 и при условии, что внешнее поле ориентиро-
вано вдоль оси oz: H0 = (0, 0, H0). Видно, что с ростом
0,04
N точки графика стремятся уложиться на единую
кривую. При наиболее возможном удалении от тор-
цов (около точки z = 0) значение нормальной состав-
ляющей находится вблизи нуля, что согласуется с
очевидным результатом решения задачи о бесконеч-
0
ном магнитном цилиндре в постоянном продольном
его оси внешнем поле. На рис. 3 приведены графики
Вариант 1
H
n
(
z
)
для различных значений длины цилиндра, и с
Вариант 2
ее увеличением видна тенденция к снижению значе-
-0,04
Вариант 3—9
Рис. 2. Зависимость нормальной составляющей результирую-
щего поля на поверхности S1 от координаты z при различных
z
значениях плотности разбиений; h1 = h2 = 1,5; ρ0 = 0,5; μ = 99;
внешнее поле H0 = (0, 0, H0).
–1
0
1
Дефектоскопия
№ 10
2019
30
В.В. Дякин, О.В. Кудряшова, В.Я. Раевский
Рис. 3. Зависимость нормальной составляющей результирующе-
го поля на поверхности S1 от координаты z при различных значе-
ниях h1 = h2; плотность разбиения по варианту 8; ρ0 = 0,5; μ = 99;
внешнее поле H0 = (0, 0, H0).
0,04
ния функции вблизи z = 0. Исходя из этого можно
судить, что при таком внешнем поле цилиндр, про-
дольный размер которого хотя бы в два раза превы-
0
шает его диаметр, при измерениях посередине вполне
можно рассматривать как неограниченный.
Рис. 4 иллюстрирует результаты расчета значения
нормальной составляющей результирующего поля в
сечении z = 0 при условии, что внешнее поле имеет
-0,04
вид H0 = (H0, 0, 0). Согласно точному решению соот-
z
ветствующей задачи о бесконечном цилиндре,
2
-2
-1
0
1
2
 =
ϕ)
cosϕ[1]. Из рисунка видно, что точки
n
µ+1
укладываются на эту кривую, разве что вблизи экс-
тремумов заметно небольшое отличие, которое можно объяснить влиянием торцов цилиндра:
на рис. 5 в увеличенном масштабе приводятся графики тех же зависимостей, построенные для
различных значений величин h1 и h2 (h1 = h2), и видно, что с их увеличением расчетная кривая при-
ближается к теоретической.
0,022
0,02
0,02
0,01
0
0,018
-0,01
Бесконечный цилиндр
Вариант 1
0,016
h1 = h2 = 3,0
Вариант 2
h1 = h2 = 2,0
h1 = h2 = 1,5
-0,02
Вариант 3—9
h1 = h2 = 1,0
φ
0,014
φ
0
2
4
6
0
0,2
0,4
0,6
0,8
Рис.
4. Зависимость нормальной составляющей
Рис. 5. Зависимость нормальной составляющей резуль-
результирующего поля на поверхности S1 от координа-
тирующего поля на поверхности S1 от координаты φ
ты φ в плоскости z = 0 при различных значениях плот-
в плоскости z = 0 при плотности разбиений по варианту
ности разбиений; h1 = h2 = 1,5; ρ0 = 0,5; μ = 99; внешнее
8 и при различных значениях h1 = h2; ρ0 = 0,5; μ = 99;
поле H0 = (H0, 0, 0).
внешнее поле H0 = (H0, 0, 0).
3. ВЫЧИСЛЕНИЕ РЕЗУЛЬТИРУЮЩЕГО ПОЛЯ
Следующий этап численных расчетов — получение значения результирующего поля. Для этого
на основе (1) с учетом принятых правил дискретизации записываем выражение
3
N
i
µ-1
dS
0
im
3
H r)=
H r)
H
(
r
)
,
r
\
S,
(12)
∑∑
n
i
im
4π
r-r
i=1
m=1
S
im
исходя из которого, после вычисления производных, выписываем окончательные формулы для
расчета компонент результирующего в цилиндрических координатах:
Дефектоскопия
№ 10
2019
К расчету поля конечного магнитного цилиндра
31
3
N
i
0
µ-1
H
τ
(r)
=
H
τ
(r)
∑∑
H
n
(
r
im
)
J
τ;im
(
r
)
,
(13)
i
4π
i
=1
m=1
где под символом τ понимается ρ, φ или z.
Выражения для Jτ;im(r) получены снятием интегрирования по одной из переменных:
1m ϕ
ϕ
J
ρ;1m
(r)
= -ρ
0
dϕ′H
(
ϕ′
)(
ρ-ρ
0
cos
(
ϕ-ϕ′
))
;
1m ϕ
ϕ
−δ
1m ϕ
ϕ
2
J
ϕ;1m
(r)
= -ρ
0
dϕ′
H
(
ϕ′
)
sin
(
ϕ-ϕ′
)
;
1m ϕ
ϕ
−δ
t=z
1m
z
z
ϕ
1m
ϕ
1
J
z;1m
(r)
= -ρ
0
d
ϕ′
,
2
2
2
ρ
2
ρρ
cos
(
ϕ-ϕ′
)
+
t
ϕ
1m
−δ
ϕ
0
0
t=z
1m
z−δ
z
где
t=z
1m
z
z
t
H
(
ϕ′
)
:=
2
2
2
2
2
ρ
2ρρ
cos
(
ϕ-ϕ′
ρ
2
ρρ
cos
(
ϕ-ϕ′
)
+
t
(
0
0
))
0
0
t=z
1m
z
−δ
z
Для i = 2, 3 при η2 = -h1 и η3 = h2 имеем
ϕ
im
ϕ
ρ′=ρ
im
ρ
J
ρ;
im
(r)
=
dϕ′F
(
ρ′
,
ϕ′ ρ,ϕ,z
,
η
i
)
,
ρ
im
ρ
ϕ
im
−δ
ϕ
где функция F определяется по формуле (10). Далее,
ϕ
1m
ϕ
ρ′=ρ
im
ρ
J
ϕ;im
(r)
=
d
ϕ′
sin
(
ϕ-ϕ'
)
G
(
ρ′
,
ϕ′ ρ,ϕ,z,
η
i
)
,
ρ′=ρ
im
ρ
ϕ
1m ϕ
−δ
где функция определяется по формуле (11). И, наконец,
ρ′=ρ
im
ρ
ϕ
im
ϕ
(
ρ′-ρ
cos
(
ϕ-ϕ
))
d
ϕ
J
(r)
:=-
(
z
−η
)
z;im
i
2
2
2
2
2
2
ϕ
ρ
sin
(
ϕ-ϕ′
)
+
(
z
−η
)
ρ
+ ρ′
+
(
z
−η
)
2
ρρ′
cos
(
ϕ-ϕ′
)
im
ϕ
(
i
)
i
im ρ
′ ρ
−δ
Для получения итогового значения приближенной величины результирующего поля по форму-
ле (13) каждый из данных определенных интегралов заменялся на произведение значения подын-
тегральной функции в точке φim с длиной 2δφ отрезка интегрирования.
На основе результатов численных экспериментов построены различные зависимости.
 ρϕ
На рис. 6а, б представлены графики (ρ,
ρ
ϕ,z)
и
z
z
)
соответственно, полученные в
результате расчетов при ρ = const и φ = 0 с внешним полем, направленным вдоль оси oz. Что каса-
ется функции (ρ,
,z),
ϕ
ϕ
то при таком внешнем поле она обращается в ноль, что подтверждают
и расчеты. Точки экстремума-компоненты и точки перегиба z-компоненты вполне четко указывают
на продольный размер цилиндра.
Если точки, в которых вычисляются значения компонент результирующего поля, располагают-
ся вблизи цилиндра, то вблизи экстремумов проявляется влияние выбранной плотности разбие-
ний, что наглядно иллюстрирует рис. 7а: расположенная выше группа кривых
полу-
ченная для вариантов разбиений с пятого по восьмой, существеннее отражает упомянутое влияние
по сравнению с аналогичной группой кривых
По мере отдаления от цилиндра абсо-
лютные значения компонент поля закономерно уменьшаются (рис. 7б и 8).
Дефектоскопия
№ 10
2019
32
В.В. Дякин, О.В. Кудряшова, В.Я. Раевский
а
б
3
h1 = 2,0
h1 = 2,0
2
h1 = 1,5
h1 = 1,5
h1 = 1,0
1
2
h1 = 1,0
0
1
-1
0
-2
z
z
–3
-2
-1
0
1
2
3
-3
-2
-1
0
1
2
3
Рис. 6. а — зависимость ρ-компоненты результирующего поля от координаты z при ρ = 0,6 и при различных значениях
h1 = h2; ρ0 = 0,5; μ = 99; внешнее поле H0 = (0, 0, H0); б — зависимость z-компоненты результирующего поля от коорди-
наты z при ρ = 0,6 и при различных значениях h
= h2; ρ0 = 0,5; μ = 99; внешнее поле H0 = (0, 0, H0).
1
а
б
3,5
8
7
ρ = 0,51
3
3
6
ρ = 0,60
ρ = 0,80
2,5
5
ρ = 1,00
2
2
1,5
1
ρ = 0,51
8
1
ρ = 0,60
5
0
z
z
0,5
1,1
1,2
1,3
1,4
1,5
1,6
1,7
0
0,4
0,8
1,2
1,6
2
Рис. 7. а — зависимость ρ-компоненты результирующего поля от координаты z при различных значениях ρ и при h1 = h2=
= 1,5 вблизи z = 1,5 для вариантов разбиений от 5 до 8; ρ
= 0,5; μ = 99; внешнее поле H0 = (0, 0, H0); б — зависимость
0
ρ-компоненты результирующего поля от координаты z при различных значениях ρ и при h1 = h2 = 1,5 вблизи z = 1,5 для
варианта разбиения 8; ρ0 = 0,5; μ = 99; внешнее поле H0 = (0, 0, H0).
Задача о бесконечном круговом цилиндре с постоянной магнитной проницаемостью, помещен-
ном в перпендикулярное его оси однородное внешнее поле, имеет простое решение в элементар-
ных функциях. В цилиндрических координатах [1]:
2
2
ρ
ρ
0
0
0
0
H
(
ρ,ϕ
)
=
H
1
cosϕ;
H
(
ρ,ϕ
)
=
H
1
sinϕ;
H
=
0.
(14)
ρ
2
ϕ
2
z
ρ
ρ
В декартовых координатах [1]:
2
2
0
2
x
y
0
2
xy
H
x,
y
=
H
1
+ χρ
;
H
x,
y
=
2H
χρ
x
(
)
0
2
y
(
)
0
2
(15)
2
2
2
2
(
x
+
y
)
(
x
+
y
)
Дефектоскопия
№ 10
2019
К расчету поля конечного магнитного цилиндра
33
5
4
ρ = 0,51
ρ = 0,60
3
ρ = 0,80
ρ = 1,00
2
1
0
-1
z
-4
-3
-2
-1
0
1
2
3
4
Рис. 8. Зависимость z-компоненты результирующего поля от координаты z при различных значениях ρ при h1 = h2 = 1,5
для варианта разбиения 8; ρ0 = 0,5; μ = 99; внешнее поле H0 = (0,0, H0).
Как и следовало ожидать, с увеличением продольного размера цилиндра кривые для компонент
результирующего поля, построенные в равноудаленной от обоих концов плоскости, имеют тенден-
цию сходиться к аналитической кривой. Рис. 9а и б отражают эту закономерность: здесь на фоне
графиков функций
H
(
x,0)
и
H
(
x,0)
даны расчетные кривые, соответствующие различным
x
y
значениям h1 = h2.
а
б
0,4
1
0,2
0,8
h1 = h2 = 1,5
h1 = h2 = 2,0
0
h1 = h2 = 1,5
h1 = h2 = 2,5
h1 = h2
= 3,0
h1 = h2 = 2,0
0,6
h1 = h2 = ∞
h1 = h2 = 2,5
-0,2
h1 = h2 = 3,0
h1 = h2
= ∞
0,4
-0,4
x
x
-1,5
-1
-0,5
0
0,5
1
1,5
–1,5
-1
-0,5
0
0,5
1
1,5
Рис. 9. а — зависимость x-компоненты результирующего поля от координаты x при y = 0,6 и z = 0 при различных зна-
чениях h1 = h2 на фоне кривой для бесконечного цилиндра; ρ0 = 0,5; μ = 99; внешнее поле H0 = (H0, 0, 0); б — зависимость
y-компоненты результирующего поля от координаты x при y = 0,6 и z = 0 при различных значениях h
= h2 на фоне кривой
1
для бесконечного цилиндра; ρ0 = 0,5; μ = 99; внешнее поле H0 = (H0, 0, 0).
К обратной задаче магнитостатики следует относится с большой осторожностью, поскольку
даже в узкой постановке она может не иметь единственного решения — эта проблема достаточно
подробно раскрыта в работе [12]. Однако, рассматривая полученные кривые зависимостей компо-
нент поля от координат, видим, что они существенно отличаются по своей форме от аналогичных
кривых для шара. Таким образом, по результатам измерений вполне возможно вынести суждение
о форме магнитного тела и даже об его геометрических размерах, если к нему нет видимого досту-
па. Детальное рассмотрение этого вопроса оставим на следующую публикацию.
Другая часть работы будет связана с развитием алгоритма при усложнении формы магнетика.
Как минимум, планируется рассмотреть цилиндр внутренним с дефектом, используя при этом
результаты настоящей работы как эталон.
Дефектоскопия
№ 10
2019
34
В.В. Дякин, О.В. Кудряшова, В.Я. Раевский
Работа выполнена в рамках государственного задания МИНОБРНАУКИ России (тема
«Квант», номер госрегистрации АААА-А18-118020190095-4).
СПИСОК ЛИТЕРАТУРЫ
1. Неразрушающий контроль и диагностика / Справочник под ред. В.В. Клюева. М.: Машиностроение,
1995. 487 с.
2. Дякин В.В., Сандовский В.А. Теория и расчет накладных вихретоковых преобразователей. М.:
Наука, 1981. 136 с.
3. Дякин В.В., Умергалина О.В., Сандовский В.А. Точное решение одной задачи магнитостатики в
бицилиндрических координатах // Дефектоскопия. 1999. № 6. С. 29—35.
4. Дякин В.В., Кудряшова О.В., Раевский В.Я. Точное решение одной задачи магнитостатики в бипо-
лярных координатах (продолжение) // Дефектоскопия. 2016. № 12. С. 48—58.
5. Дякин В.В., Кудряшова О.В. Дефект в цилиндре // Дефектоскопия. 2012. № 4. С. 41—55.
6. Дякин В.В., Кудряшова О.В. Дефект в трубе // Дефектоскопия. 2012. № 10. С. 3—17.
7. Philippe Ciarlet. The Finit Element Method for Elliptic Problems. North-Holland Publishing Company,
Amsterdam — New York — Oxford, 1978. 512 p.
8. Хижняк Н.А. Интегральные уравнения макроскопической электродинамики. Киев: Наукова
думка, 1986. 278 с.
9. Михлин С.Г. Курс математической физики. М.: Наука, 1968. 575 с.
10. Градштейн И.С., Рыжик И.М. Таблицы интегралов, сумм, рядов и произведений. М.: ГИФМЛ,
1962. 1100 с.
11. George E. Forsythe, Michael A. Malcolm, Cleve B. Mouler. Computer Methods for Mathematical
Computations. PRENTICE-HALL, INC, ENGLWOOD CLIFFS, N.J. 07632, 1977. 279 p.
12. Dyakin V. V., Kydryashova O.V., Raevskii V.Y. On the Well-Posedness of the Direct and Inverse Problem
of Magnetostatics. Part 2 // Russian Journal of Nondestructive testing. 2018. V. 54. No. 10. P. 687—697.
[Дякин В.В., Кудряшова О.В., Раевский В.Я. К вопросу корректности прямой и обратной задачи магни-
тостатики. Часть 2 // Дефектоскопия. 2018. № 10. С. 15—24.]
Дефектоскопия
№ 10
2019