How to model permanent magnets in MATLAB pde?

I am trying to model the magnetic flux density distribution of permanent magnets in MATLAB pde, but I have faced a serious problem. MATLAB pde can consider the elliptic pde of the form: -1/Mu0*del2(A)=J, where J is the free current density in the considered area. However, here we have a couple of permanent magnets (PMs, N and S) with constant magnetization vectors (M1=Br/Mu0 and M2=-Br/Mu0) inside it, where Br is the residual magnetic flux density of PMs. So the area of PMs can be modelled with a periodic pulse wise magnetization vector M. A conventional answer is to model the PMs as some equivalent bounded current densities Jb=curl(M), which when M=My(x)j, we have Jb=dMy/dx. However, because of the pulse wise nature of M, the differentiation gives infinite values for Jb at the PM corners. How can I solve this problem? I will greatly appreciate someone who can help me.

[Hans DeSmidt]

Hi, I am facing the same issue.

[Volker Ziemann]

You can model permanent magnets through surface currents. For the theory you can check out section 4.5 in V. Ziemann, Hands on Accelerator Physics Using Matlab, CRC Press. The code below illustrates how to model a permanent magnet cube with the PDE toolbox (works on R2022b). Optionally you can replace the surface currents by current sheets. Their field can be calculated using Biot-Savart and this is described in https://arxiv.org/abs/2106.04153. I used that for teaching and the code is available from https://github.com/volkziem/PermanentMagnetDesign.

% PM_cube.m, V. Ziemann, 231206
clear all; close all
model=createpde(1);
%..........................geometry
Br=1;     % remanent field in Tesla
dw=0.01;  % width of current sheet
x0=0.5;   % half width of cube
y0=x0;    % half height of cube
CS1=[2;4; x0-dw/2;x0+dw/2;x0+dw/2;x0-dw/2;-y0;-y0;y0;y0]; 
CS2=[2;4; -x0-dw/2;-x0+dw/2;-x0+dw/2;-x0-dw/2;-y0;-y0;y0;y0]; 
ws=5;     % simulation volume, box
World=[2;4; -ws;ws;ws;-ws;-ws;-ws;ws;ws]; 
gd=[World,CS1,CS2];             % assemble geometry
ns=char('World','CS1','CS2')';  % names of the regions                 
sf='World+CS1+CS2'; 
[dl,bt]=decsg(gd,sf,ns);
geometryFromEdges(model,dl);
figure(1); clf;
pdegplot(model,'EdgeLabels','on'); axis off; axis square;
figure(2); clf; 
pdegplot(model,'SubDomainLabels','on'); axis square;
%return   % check the assignment of edges and faces
generateMesh(model,'Hmax',0.1);
figure(3); clf; pdemesh(model); axis square;
applyBoundaryCondition(model,'Edge',[1,2,4,5],'u',0);
specifyCoefficients(model,'m',0,'d',0,'c',1,'a',0,'f',-Br/dw,'Face',3);
specifyCoefficients(model,'m',0,'d',0,'c',1,'a',0,'f',Br/dw,'Face',2);
specifyCoefficients(model,'m',0,'d',0,'c',1,'a',0,'f',0,'Face',1);
result=solvepde(model);
figure(4); clf; pdeplot(model,'xydata',result.NodalSolution,'contour','on'); 
set(gca,'FontSize',16)
hold on; pdegplot(model)
By=result.XGradients; Bx=-result.YGradients; 
figure(5); clf; pdeplot(model,'xydata',result.NodalSolution,'xystyle','off', ...
    'flowdata',[Bx,By],'contour','on','colorbar','off'); axis equal;
set(gca,'FontSize',16)
y=-5*x0:dw/2:5*x0; x=0*ones(1,length(y));  % along vertical line
[By,Bx]=evaluateGradient(result,x,y);  Bx=-Bx;
B=hypot(Bx,By);  % absolute value of B
figure(6); plot(y,B,'k','LineWidth',2);
xlabel('y [m]'); ylabel('B [T]')
set(gca,'FontSize',16)