Расчёт полей тороидально намагниченного цилиндра

[liquidcrystalosmos]

Представляю на суд публики свой расчёт полей тороидально намагниченного цилиндра

Расчёт выполнен в среде с открытым исходным кодом SageMath-9.1

Данная среда позволяет сделать экспорт расчёта в виде веб-страницы, чем я и воспользовался

На сегодняшний день я выполнил два варианта расчёта:

первый вариант - модель постоянной намагниченности
http://liquidcrystalosmos.narod.ru/field_of_deyna_cylinder.htm

второй вариант - модель, в которой намагниченность обратно пропорциональна радиусу
http://liquidcrystalosmos.narod.ru/field_of_deyna_cylinder_with_magnetization_inverse_to_radius.htm

В основу расчёта я взял понятие векторного потенциала. Пространственные производные которого (ротор и дивергенция) суть магнитные поля.

В чем состоит вопрос. Удивительным для меня образом я получил отличную от нуля дивергенцию векторного потенциала, хотя, например, в курсе электродинамики Тамма доказывается, что для замкнутых постоянных токов таковая должна быть равна нулю.

В чем состоит ошибка, если таковая имеется?

[meandr]

Эти эксперименты Николаева и Дейны обсуждались на Сайтехе
http://www.sciteclibrary.ru/cgi-bin/yabb2/YaBB.pl?num=1510731773
Знатоки пришли к выводу, что весь эффект - из-за нарушений ТОЧНОЙ симметричной тороидальности поля
- либо стальные кольца сдвинуты относительно оси = намагничивающего провода с током,
- либо в самом материале колец неоднородности, приводящие к неоднородности поля
и т.п.

Я со своей стороны замечу, что разногласия получаются при использовании "токовой" модели самого намагниченного тела - протекающими по его поверхности "токами Ампера". В этой модели вроде бы должно получаться взаимодействие даже идеально симметричных тороидальных магнетиков="тороидальных обмоток с током" - но на самом деле для идеального случая с магнетиками такого взаимодействия не должно быть.

В основу расчёта я взял понятие векторного потенциала. Пространственные производные которого (ротор и дивергенция) суть магнитные поля.
В чем состоит вопрос. Удивительным для меня образом я получил отличную от нуля дивергенцию векторного потенциала, хотя, например, в курсе электродинамики Тамма доказывается, что для замкнутых постоянных токов таковая должна быть равна нулю.

В чем состоит ошибка, если таковая имеется?

Насчет того что индукция B=rotA я в курсе.
А вот что такое магнитное у ТС подразумевает divA ?

Обычно пишут либо в калибровке Кулона divA=0 всегда и везде, либо в калибровке Лоренца
\[ c^2divA=-d\varphi/dt \]
тогда для статического случая тоже divA=0 как у Тамма так и у других Учителей.

Не влазя в подробности Ваших весьма объемных вычислений, не могу сказать, где именно Вы ошиблись, скорее всего где-то написав "следуя Тамму" ВЫ на самом деле ему не последовали (или сам Тамм туда не пошел бы).
Возьмем случай, когда намагниченность М в тороиде (и индукция В вместе с ней) убывают обратно пропорционально радиусу (как в поле вокруг длинного прямого проводника с током).
Можете указать, в каком конкретно месте такого тороида у Вас вмп А зарождается (divA>0) и где исчезает (divA<0) ?

[liquidcrystalosmos]

Возьмем случай, когда намагниченность М в тороиде (и индукция В вместе с ней) убывают обратно пропорционально радиусу (как в поле вокруг длинного прямого проводника с током).
Можете указать, в каком конкретно месте такого тороида у Вас вмп А зарождается (divA>0) и где исчезает (divA<0) ?

Полученные мною графики дивергенции векторного потенциала А (взятые со знаком минус) выглядят вот так


Интересен также нормированный график ротора Н

[liquidcrystalosmos]

Геометрия задачи


Рассмотрим случай в котором намагниченность обратно пропорциональна радиусу. Этот случай проще потому что объёмные токи равны нулю.

Цилиндр расположен у меня горизонтально, как на рисунке, поэтому торцы у меня правый и левый.

расчёты я делаю в сгс

J - ток по прямолинейному медному проводу в опыте Дейны

напряжённость магнитного поля прямого тока H = 2 * J / (c * b)
   источник: Савельев, т.2  стр.423

вектор намагниченности равен I = kappa * H = kappa * 2 * J / (c * r)

поверхностный ток у Тамма \(js = c \cdot [I \times n]\) (параграф 61.)

Код: [Выделить]

c = var("c")                   # скорость света
kappa = var("kappa")           # Магнитная восприимчивость
J = var("J")                   # ток
js(J, kappa, r) = 2*J*kappa/r  # плотность поверхностного тока на внутренней и внешней цилиндрической поверхностях цилиндра
                               # js(rj1) = - 2*J*kappa/rji
                               # js(rj2) = + 2*J*kappa/rj2
jt(J, kappa, r) = 2*J*kappa/r  # плотность поверхностного тока на торцах цилиндра
                               # jt(rj) = + 2*J*kappa/rj # zj1
                               # jt(rj) = - 2*J*kappa/rj # zj2

на внешней цилиндрической поверхности направление тока положительно, на внутренней - отрицательно.
на левой торцевой поверхности направление тока положительно, на - правой отрицательно.

Код: [Выделить]

z_j1 = var("z_j1") # левый торец цилиндра
z_j2 = var("z_j2") # правый торец цилиндра

r_j1 = var("r_j1") # радиус внутренней поверхности (сверления) цилиндра
r_j2 = var("r_j2") # радиус внешней цилиндрической поверхности цилиндра


[liquidcrystalosmos]

Интеграл по \({\varphi}_{j}\) векторного потенциала, создаваемого элементом площади \({{r}_{j}}\,d{{\varphi}_{j}}d{z_j}\) цилиндрической поверхности цилиндра

\[A_s=\frac{1}{c}\int\limits_{{{\varphi}_{j}}}{\frac{j_s \left( {{r}_{j}} \right){{r}_{j}}}{R}\ }d{{\varphi}_{j}}\]

Интеграл по \({\varphi}_{j}\) векторного потенциала, создаваемого элементом площади \({{r}_{j}}\,d{{\varphi}_{j}}d{r_j}\) торцевой поверхности цилиндра

\[A_t=\frac{1}{c}\int\limits_{{{\varphi}_{j}}}{\frac{j_t \left( {{r}_{j}} \right){{r}_{j}}}{R}\ }d{{\varphi }_{j}}\]

интеграл \(IR_{\varphi}=\int\limits_{0}^{2\pi}{\frac{1}{R}}d{\varphi_j}\) выражается через полный эллиптический интеграл первого рода следующим образом:

Код: [Выделить]

rja2   = lambda rj, ra, zj, za : (rj-ra)^2+(zj-za)^2
module = lambda rj, ra, zj, za : - 4*rj*ra / rja2(rj, ra, zj, za)
IRphi  = lambda rj, ra, zj, za : 4*elliptic_kc(module(rj, ra, zj, za)) / sqrt(rja2(rj, ra, zj, za))

\[IR_{\varphi} = \frac{4 \, K\left(-\frac{4 \, r_{a} r_{j}}{{\left(r_{a} - r_{j}\right)}^{2} + {\left(z_{a} - z_{j}\right)}^{2}}\right)}{\sqrt{{\left(r_{a} - r_{j}\right)}^{2} + {\left(z_{a} - z_{j}\right)}^{2}}}\]

Код: [Выделить]

As = lambda J, c, kappa, rj, ra, zj, za : IRphi(rj, ra, zj, za) * js(J, kappa, rj)*rj/c
At = lambda J, c, kappa, rj, ra, zj, za : IRphi(rj, ra, zj, za) * jt(J, kappa, rj)*rj/c

\[As = \frac{8 \, J \kappa K\left(-\frac{4 \, r_{a} r_{j}}{{\left(r_{a} - r_{j}\right)}^{2} + {\left(z_{a} - z_{j}\right)}^{2}}\right)}{\sqrt{{\left(r_{a} - r_{j}\right)}^{2} + {\left(z_{a} - z_{j}\right)}^{2}} c}\]

\[At = \frac{8 \, J \kappa K\left(-\frac{4 \, r_{a} r_{j}}{{\left(r_{a} - r_{j}\right)}^{2} + {\left(z_{a} - z_{j}\right)}^{2}}\right)}{\sqrt{{\left(r_{a} - r_{j}\right)}^{2} + {\left(z_{a} - z_{j}\right)}^{2}} c}\]

[liquidcrystalosmos]

Рассчитаем производную векторного потенциала, создаваемого торцевым поверхностным током, по \(z_a\) координате точки наблюдения \(\frac{\partial}{\partial z_a} A_t = \frac{\partial}{\partial z_a}\int\limits_{{{\varphi}_{j}}}{\frac{j_t \left( {{r}_{j}} \right){{r}_{j}}}{R}\ }d{{\varphi }_{j}}\)

Код: [Выделить]

At_diff_za_ = lambda J, c, kappa, rj, ra, zj, za : diff(At(J, c, kappa, rj, ra, zj, za), za)

Учитывая положительное направление поверхностного тока в левом торце и отрицательное направление поверхностного тока в правом торце составим сумму полученной производной для обоих торцов

Код: [Выделить]

At2_diff_za_ = lambda J, c, kappa, rj, ra, zj1, zj2, za : At_diff_za_(J, c, kappa, rj, ra, zj1, za) - At_diff_za_(J, c, kappa, rj, ra, zj2, za)

Рассчитаем производную векторного потенциала, создаваемого торцевым поверхностным током, по $r_a$ координате точки наблюдения $\frac{1}{r_a}\frac{\partial}{\partial r_a}\left(r_a\,A_t\right) = \frac{1}{r_a}\frac{\partial}{\partial r_a}\left(r_a\,\int\limits_{{{\varphi}_{j}}}{\frac{j_t \left( {{r}_{j}} \right){{r}_{j}}}{R}\ }d{{\varphi }_{j}}\right)$

Код: [Выделить]

At_ra_diff_ra_div_ra_ = lambda J, c, kappa, rj, ra, zj, za : (ra*At(J, c, kappa, rj, ra, zj, za)).diff(ra)/ra

Учитывая положительное направление поверхностного тока в левом торце и отрицательное направление поверхностного тока в правом торце составим сумму полученной производной для обоих торцов

Код: [Выделить]

At2_ra_diff_ra_div_ra_ = lambda J, c, kappa, rj, ra, zj1, zj2, za : At_ra_diff_ra_div_ra_(J, c, kappa, rj, ra, zj1, za) - At_ra_diff_ra_div_ra_(J, c, kappa, rj, ra, zj2, za)

Рассчитаем производную векторного потенциала, создаваемого поверхностным током цилиндрической поверхности, по $r_a$ координате точки наблюдения $\frac{\partial}{\partial r_a} A_s = \frac{\partial}{\partial r_a}\int\limits_{{{\varphi}_{j}}}{\frac{j_s \left( {{r}_{j}} \right){{r}_{j}}}{R}\ }d{{\varphi }_{j}}$

Код: [Выделить]

As_diff_ra_ = lambda J, c, kappa, rj, ra, zj, za : As(J, c, kappa, rj, ra, zj, za).diff(ra)

Учитывая отрицательное направление поверхностного тока на внутренней цилиндрической поверхности (в сверлении) и положительное направление поверхностного тока на внешней цилиндрической поверхности составим сумму полученной производной для обоих цилиндрических поверхностей

Код: [Выделить]

As2_diff_ra_ = lambda J, c, kappa, rj1, rj2, ra, zj, za : - As_diff_ra_(J, c, kappa, rj1, ra, zj, za) + As_diff_ra_(J, c, kappa, rj2, ra, zj, za)

Рассчитаем производную векторного потенциала, создаваемого поверхностным током цилиндрической поверхности, по $z_a$ координате точки наблюдения $\frac{\partial}{\partial z_a} A_s = \frac{\partial}{\partial z_a}\int\limits_{{{\varphi}_{j}}}{\frac{j_s \left( {{r}_{j}} \right){{r}_{j}}}{R}\ }d{{\varphi }_{j}}$

Код: [Выделить]

As_diff_za_ = lambda J, c, kappa, rj, ra, zj, za : As(J, c, kappa, rj, ra, zj, za).diff(za)

Учитывая отрицательное направление поверхностного тока на внутренней цилиндрической поверхности (в сверлении) и положительное направление поверхностного тока на внешней цилиндрической поверхности составим сумму полученной производной для обоих цилиндрических поверхностей

Код: [Выделить]

As2_diff_za_ = lambda J, c, kappa, rj1, rj2, ra, zj, za : - As_diff_za_(J, c, kappa, rj1, ra, zj, za) + As_diff_za_(J, c, kappa, rj2, ra, zj, za)


[liquidcrystalosmos]

Дальнейшие пояснения по методике расчёта

1) Замечаю, что ввиду цилиндрической симметрии задачи $\frac{\partial A}{\partial \varphi_a} = 0$ и мне не потребуется дифференцирование векторного потенциала по $\varphi_a$ координате точек наблюдения, выражаю интеграл от обратного радиуса по $\varphi_j$ координате точек истока $IR_{\varphi}=\int\limits_{0}^{2\pi}{\frac{1}{R}}d{\varphi_j}$ через полный эллиптический интеграл первого рода

2) умножая полученное выше выражение на плотность поверхностного тока получаю выражение для векторного потенциала создаваемого кольцом на поверхности бесконечно малой ширины.

В обозначениях моего расчёта:

векторный потенциал бесконечно тонкого кольца торцевой поверхности At
векторный потенциал бесконечно тонкого кольца цилиндрической образующей поверхности As

3) дифференцирую векторные потенциалы колец по координатам точек наблюдения и получаю аналитическое выражение частных производных векторного потенциала поверхностного тока колец по координтам точек наблюдения

At_diff_za
At_ra_diff_ra_div_ra
As_diff_ra
As_diff_za

4) Суммирую эти производные (учитывая знаки направления токов) по смежным парам колец, получаю

At2_diff_za
At2_ra_diff_ra_div_ra
As2_diff_ra
As2_diff_za

5) Подставляю численные значения которые можно подставить на данном этапе (координаты торцов и радиусов, ток J, магнная восприимчивость и т.д)

Дальнейшие действия произвожу внутри функций по вычислению магнитного поля Н и дивергенции векторного потенциала.

Эти функции на вход принимают координаты точек наблюдения. На выходе выдают поле.

Внутри функции:

6) подставляю координаты точек наблюдения.

7) Произвожу численное интегрирование полученных производных векторного потенциала колец по координатам истока в соответствии с размерами цилиндра, превращая бесконечно тонкие кольца истока в торцевые и образующие поверхности истока.

8) Суммируя полученные интегралы, нахожу поле Н и дивергенцию А в точке наблюдения

H_phi = At_diff_za_num_int - As_diff_ra_num_int
H_scalar = - div A = - At_ra_diff_ra_div_ra_num_int = - As_diff_za_num_int

[liquidcrystalosmos]

Ошибка в моих расчётах найдена http://www.sciteclibrary.ru/cgi-bin/yabb2/YaBB.pl?num=1614926067/20#20