Вычисление магнитной индукции постоянного цилиндр. магнита

[GAA]

Используется система СИ.
Пусть (помещённый в "вакуум") магнит имеет форму цилиндра $x^2+y^2 \le R^2$, $-h \le z \le 0$. Вектор (остаточной) намагниченности не зависит от координат (внутри цилиндра) и направлен вдоль оси $Oz$. Внешнее поле отсутствует. Тогда скалярный потенциал $ \mathbf {H}(q) = -\nabla \varphi (q)$ можно представить в виде

$$\varphi (q) = \frac{1}{4\pi}\int_V \mathbf{M}\nabla_r \frac{1}{| q- r|} dv,$$

где $V$ цилиндр (по которому выполняется интегрирование), $q = (x_0, y_0, z_0)$, $r = (x, y, z)$.
Т.к. $\mathbf{M}$ направлен вдоль оси $Oz$ и не зависит от координат, выражение для $\varphi$ можно записать в виде

$\varphi (q) = \frac{M}{4\pi}\int_V \partial_z \frac{1}{| q- r|} dv$ или $\varphi (q) = \frac{M}{4\pi}\int_V \frac{z_0-z}{| q- r|^3} dv.$

Здесь $M$ — модуль (длина) $\mathbf{M}$.

Интегрирование по переменной $z$ приводит к двойному интегралу

$$\varphi (q) = \frac{M}{4\pi}\left\{ \int\limits_{x^2+y^2 \le R^2} \frac {dxdy}{\sqrt {(x_0-x)^2+(y_0-y)^2+z_0^2} } - \int\limits_{x^2+y^2 \le R^2} \frac {dxdy}{\sqrt {(x_0-x)^2+(y_0-y)^2+(z_0+h)^2} } \right\}.$$

Сводя двойной интеграл к повторному и интегрируя внутренний интеграл по $y$, получим интеграл по $x$ от длинного логарифма. После вычисления градиента выражение становится довольно громоздким.

Рассмотрим значения $\mathbf H$ на оси $Oz$. На этой оси координаты $H_x$ и $H_y$ равны нулю.
При $z_0 \ne 0$, $z_0 \ne -h$

$$H_z (z_0)= \frac M 2 \left\{ z_0 \left( \frac 1 {|z_0|} - \frac 1 {\sqrt{R^2 +z_0^2}}\right) -(z_0+h) \left(\frac 1 {|z_0+h|} - \frac 1 {\sqrt {R^2 +(z_0+h)^2}} \right) \right\}.$$

$H_z(+0) = H_z(-h-0) = \frac {Mh}{2\sqrt{R^2+h^2}}$, $H_z(-0) = H_z(-h+0) = -M +\frac {Mh}{2\sqrt{R^2+h^2}}$.
Зависимость $H_z$ от $z$ для $M=1$, $h=1$ и $R=1$ приведена на рис.

Вложение:

Pic1.PNG

А на следующем рисунке приведена зависимость $B_z/\mu_0 = H_z + M_z$ ($M_z = M$).
(Как и должно быть $B_z$ «непрерывна».)

Вложение:

Pic2.PNG

Разумный ли способ вычисления $\mathbf{B}$ выбран? Имеется ли способ более быстрых вычислений (не выполнять численный расчёт однократного интеграла для случая $x_0 \ne 0$, $y_0 \ne 0$)? В каких книгах можно посмотреть? ($\mathbf{B}$ планируется использовать для дальнейших вычислений. Вычислять значения $\mathbf{B}$ необходимо многократно в разных точках верхнего полупространства ($z>0$).)

[amon]

GAA в сообщении #1569996 писал(а):

$$\varphi (q) = \frac{1}{4\pi}\int_V \mathbf{M}\nabla_r \frac{1}{| q- r|} dv,$$

Это то, что Вы делаете, но чуть другими словами. Известно, что для однородно намагниченного тела скалярный магнитный потенциал $\psi$ равен

$$\psi=-\mathbf{M}\cdot\nabla\varphi,$$

где $\varphi$ - электростатический потенциал того же тела с единичной плотностью заряда. То есть если мы умеем решать электростатическую задачу, то магнитостатический ответ пишется сразу. Правда, поле равномерно заряженного цилиндра мне не известно, и, боюсь, ни через что человеческое не выражается.

[GAA]

Для равномерно электрически заряженного с поверхностной плотностью $\sigma$ круга $x^2 + y^2 \le R^2$, $z=0$ в книге [1] приводится (СГС) выражение потенциала в цилиндрической системе координат

$$V =2\pi R\sigma \int\limits_0^{\infty} J_1(R l) J_0 (r l) \exp(-z_0 l) \frac {dl} l, \qquad(z_0>0)$$

где $J_m$ — функции Бесселя первого рода.

Используя этот подход (заменяя электрически заряженный круг на магнитно заряженный и учтывая два круга: с $z=0$ и $z=-h$), в системе СИ можно записать выражение для скалярного магнитного потенциала

$$\psi = \frac {MR}2 \int_0^{\infty} J_1(R l) J_0(r l) \exp(-z_0 l) (1-\exp(-hl))\frac {dl}l \qquad (z_0 >0).$$

Дифференцируя под знаком интеграла (и учитывая при вычислении производной по $r$, что $J_0’(u) = -J_1(u)$), получим

$$H_r = - \frac {\partial \psi}{\partial r} =\frac {MR} 2 \int\limits_0^{\infty} J_1(R l) J_1(r l) \exp(-z_0 l)(1-\exp(-hl)) dl;$$

$$H_z = - \frac {\partial \psi}{\partial z_0} =\frac {MR} 2\int\limits_0^{\infty} J_1(R l) J_0(r l) \exp(-z_0 l)(1-\exp(-hl)) dl.$$

Используя результаты [2] выражения для $H_r$ и $H_z$ могут быть выражены через полные эллиптические функции $K(k^2)$ (первого рода), $E(k^2)$ (второго рода) и $\Pi(p, k^2)$ (третьего рода).
Вводя обозначения $k^2 = \frac {4Rr}{(R+r)^2 +z_0^2}$, $I(z_0) = \frac{(1-k^2/2)K(k^2)-E(k^2)}{\sqrt{k^2}Rr},$будем иметь

$$H_r = \frac {RM}{\pi}(I(z_0)- I(z_0+h)), \qquad (z_0 >0).$$

Вводя обозначения $\sin^2 \beta = \frac {z_0^2}{(R-r)^2+z_0^2}$, $p=\frac {k^2}{1-(1-k^2)\sin^2\beta}$, $\Lambda = \sqrt{1-p}\sqrt{1-k^2/p} \Pi(p, k^2)$, будем иметь

$$I(z_0) = \left\{ \begin{array} 1 1 - \frac {kz_0}{2\pi\sqrt{Rr}}K(k^2)- \frac 1 {\pi}\Lambda(p, k^2), \quad r<R; \\ \frac 1 2 - \frac {kz_0}{2\pi R}K(K^2), \qquad\qquad\qquad\quad r=R; \\ -\frac {kz_0}{2\pi\sqrt{Rr}}K(k^2)+\frac 1 {\pi} \Lambda(p, k^2), \qquad r> R. \end{array}$$

$$H_z =\frac{M}{2}(I(z_0) - I(z_0+h)), \qquad (z_0 >0).$$

Ниже приведены тексты на matlab.
Файл CHx.m:

Используется синтаксис Matlab M

 function [Sx] = CHx(x0, y0, z0, M, R, h)
 fx = @(x) (x0-x)./sqrt( (sqrt(R^2-x.^2)-y0).^2 +(x0-x).^2 + z0^2)./(sqrt(R^2-x.^2)-y0 + sqrt((sqrt(R^2-x.^2)-y0).^2 +(x0-x).^2 + z0^2));
 Sx = integral(fx, -R, R);
 fx = @(x) (x0-x)./sqrt( (sqrt(R^2-x.^2)+y0).^2 +(x0-x).^2 + z0^2)./(-sqrt(R^2-x.^2)-y0 + sqrt((sqrt(R^2-x.^2)+y0).^2 +(x0-x).^2 + z0^2));
 Sx = Sx - integral(fx, -R, R);
 fx = @(x) (x0-x)./sqrt( (sqrt(R^2-x.^2)-y0).^2 +(x0-x).^2 + (z0+h)^2)./(sqrt(R^2-x.^2)-y0 + sqrt((sqrt(R^2-x.^2)-y0).^2 +(x0-x).^2 + (z0+h)^2));
 Sx = Sx - integral(fx, -R, R);
 fx = @(x) (x0-x)./sqrt( (sqrt(R^2-x.^2)+y0).^2 +(x0-x).^2 + (z0+h)^2)./(-sqrt(R^2-x.^2)-y0 + sqrt((sqrt(R^2-x.^2)+y0).^2 +(x0-x).^2 + (z0+h)^2));
 Sx = Sx + integral(fx, -R, R);
 Sx = -M/(4*pi)*Sx;
end

Файл CHz.m:

Используется синтаксис Matlab M

 function [Sz] = CHz(x0, y0, z0, M, R, h)
 fz = @(x) z0./sqrt( (sqrt(R^2-x.^2)-y0).^2 +(x0-x).^2 + z0^2)./(sqrt(R^2-x.^2)-y0 + sqrt( (sqrt(R^2-x.^2)-y0).^2 +(x0-x).^2 + z0^2));
 Sz = integral(fz, -R, R);
 fz = @(x) z0./sqrt( (sqrt(R^2-x.^2)+y0).^2 +(x0-x).^2 + z0^2)./(-sqrt(R^2-x.^2)-y0 + sqrt( (sqrt(R^2-x.^2)+y0).^2 +(x0-x).^2 + z0^2));
 Sz = Sz - integral(fz, -R, R);
 fz = @(x) (z0+h)./sqrt( (sqrt(R^2-x.^2)-y0).^2 +(x0-x).^2 + (z0+h)^2)./(sqrt(R^2-x.^2)-y0 + sqrt( (sqrt(R^2-x.^2)-y0).^2 +(x0-x).^2 + (z0+h)^2));
 Sz = Sz - integral(fz, -R, R);
 fz = @(x) (z0+h)./sqrt( (sqrt(R^2-x.^2)+y0).^2 +(x0-x).^2 + (z0+h)^2)./(-sqrt(R^2-x.^2)-y0 + sqrt( (sqrt(R^2-x.^2)+y0).^2 +(x0-x).^2 + (z0+h)^2));
 Sz = Sz + integral(fz, -R, R);
 Sz = -M/(4*pi)*Sz;
end

Файл сценария для оценки времени выполнения этих двух подходов:

Используется синтаксис Matlab M

 a=9e-6; %m, 9 micrometer
M = 20000; % A/m
Hx0 = 0.1*M;
Hz0 = 0.1*M;
h = 6e-6; %m
D = 5.4e-6; %m
R = D/2;
mu0 = 4*pi*10^(-7);
z0=1e-6;
xx=[0:a/200:3*a];
N = length(xx);
Hxp = zeros(1, N);  Hzp = zeros(1, N);
AHxp = zeros(1, N); AHzp = zeros(1, N);
tic;
for i=1:1:N
 r = xx(i);
 sqr_k = 4*R*r/((R+r)^2+z0^2);
 [K, E] = ellipke(sqr_k);
 AHxp(i) = R*M*((1-sqr_k/2)*K-E)/sqrt(sqr_k*r*R)/pi;
 sqrsinBeta = z0^2/((R-r)^2+z0^2);
 p = sqr_k/(1-(1-sqr_k)*sqrsinBeta);
 L = sqrt(1-p)*sqrt(1-sqr_k/p)*ellipticPi(p, sqr_k); 
 
 if R > r
  if r == 0
   AHzp(i) = -z0/sqrt(R^2+z0^2);   
  else    
  AHzp(i) = 1- sqrt(sqr_k)*z0/(2*pi*sqrt(R*r))*K - L/pi; 
  end
 elseif R==r
  AHzp(i) = 0.5 - sqrt(sqr_k)*z0/(2*pi*R)*K;  
 else
  AHzp(i) = -sqrt(sqr_k)*z0/(2*pi*sqrt(R*r))*K + L/pi;   
 end

 sqr_k = 4*R*r/((R+r)^2+(z0+h)^2);
 [K, E] = ellipke(sqr_k);
 AHxp(i)= AHxp(i) - R*M*((1-sqr_k/2)*K-E)/sqrt(sqr_k*r*R)/pi;
 sqrsinBeta = (z0+h)^2/((R-r)^2+(z0+h)^2);
 p = sqr_k/(1-(1-sqr_k)*sqrsinBeta);
 L = sqrt(1-p)*sqrt(1-sqr_k/p)*ellipticPi(p, sqr_k);
 if R > r
  if r==0
   AHzp(i) = M*(AHzp(i)+ (z0+h)/sqrt(R^2+(z0+h)^2))/2;      
  else    
  AHzp(i) = M*(AHzp(i)- (1- sqrt(sqr_k)*(z0+h)/(2*pi*sqrt(R*r))*K - L/pi))/2; 
  end
 elseif R==r
  AHzp(i) = M*(AHzp(i)- (0.5 - sqrt(sqr_k)*(z0+h)/(2*pi*R)*K))/2;  
 else
  AHzp(i) = M*(AHzp(i) - (-sqrt(sqr_k)*(z0+h)/(2*pi*sqrt(R*r))*K + L/pi))/2;   
 end
end
timerVal1 = toc;
tic;
for i=1:1:N
 r = xx(i);
 Hxp(i) = CHx(r, 0, z0, M, R, h);
 Hzp(i) = CHz(r, 0, z0, M, R, h);
end
timerVal2 = toc;
disp(timerVal1/timerVal2);
plot(xx, AHxp, 'b-', xx, Hxp, 'r--', xx, AHzp, 'g-', xx, Hzp, 'r-.');
xlabel('r, m'); ylabel('H_i, A/m');
legend('analytical H_x', 'Numerical H_x', 'analytical H_z', 'Numerical H_z');
title(strcat('R=', num2str(R), 'm,  M=', num2str(M), 'A/m,  z_0=', num2str(z0), 'm'));
figure;
plot(xx, AHxp-Hxp, xx, AHzp-Hzp);
legend('analytical H_x - Numerical CH_x',  'analytical H_z - Numerical H_z');
title(strcat('R=', num2str(R), 'm,  M=', num2str(M), 'A/m,  z_0=', num2str(z0), 'm'));

В R2013b время выполнения численного интегрирования приблизительно в 8 раз меньше времени выполнения кода на основе сведения к эллиптическим функциям (на очень древнем процессоре), а в R2022 приблизительно в 15 раз (на более позднем процессоре).

На рис. ниже приведены графики зависимости $H_r$ и $H_z$ от $r$ для фиксированного значения $z_0$. Абсолютная разница не более $4 \cdot 10^{-9}$.

Вложение:

HxAndHz.png

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

Ref
1. Bateman H. Partial differential equations of mathematical physics — N.Y., Dover publications, 1944; Ch VII Cylindrical co-ordinates, §7.32 The use of definite integrals involving Bessel functions, c. 409–410.
2. Eason G., Noble B., Sneddon I. N. On Certain Integrals of Lipschitz-Hankel Type Involving Products of Bessel Functions // Phil. Trans. R. Soc. Lond. A, 1955, 247.

[peregoudov]

Я это когда-то считал, может, будет полезно

[Ignatovich]

amon в сообщении #1570011 писал(а):

Это то, что Вы делаете, но чуть другими словами. Известно, что для однородно намагниченного тела скалярный магнитный потенциал $\psi$ равен

$$\psi=-\mathbf{M}\cdot\nabla\varphi,$$

где $\varphi$ - электростатический потенциал того же тела с единичной плотностью заряда.

Наверное, не единичной плотностью заряда, а с единичной поляризацией.