Hello,
I downloaded the 2D example of the magnetostatic problem at https://modules.freefem.org/modules/magnetostatic/ and I have added a few lines to compute By = dy(Az) and plot its values in 2D. According to the sign of the currents the value of By at the center should be positive, but it turns out to be negative when i run the script. Also, if I plot the current obtained from J = dx(Hy) - dy(Hx) to check for self-consistency, i see that the currents are exactly opposite with respect to the ones set at the beginning. Also, I think the signs of the weak form of the equation are correct.
Could somebody explain this behavior to me?
Thank you very much,
Adele
Shouldn’t it read Bx=-dy(Az), or By=dx(Az)?
Denis
Dear Denis,
thank you for pointing out the typo in my message.
Because B = curl(A) the correct way of writing is Bx = dy(Az).
The other component is By = - dx(Az).
I copied my script and I pasted it below. If you run it, you will se that the currents computed to check self-consistency are opposite with respect to the currents imposed as boundary conditions.
I don’t know why it comes out this way, they should be exactly the same.
Thank you,
Alberto
//Parameters
real Mu0 = 4.pi1.e-7; //Vacuum magnetic permeability
real MuC = 1.25e-6; //Copper magnetic permeability
real J0 = 500.e4; //Current density
//Mesh
int nn = 25; //Mesh quality
real BoxRadius = 2.; //“Inifite” boundary
real MagnetHeight = 1.; //Magnet height
real MagnetInternalRadius = 0.1; //Magnet internal radius
real MagnetExternalRadius = 0.3; //Magnet external radius
int BoxWall = 10;
int nnL = max(2., BoxRadiusnn/2.);
border b1(t=0., 2.pi){x=BoxRadiuscos(t); y=BoxRadiussin(t); label=BoxWall;};
border b2(t=0., 1.){x=MagnetInternalRadius+(MagnetExternalRadius-MagnetInternalRadius)t; y=-MagnetHeight/2.;};
border b3(t=0., 1.){x=MagnetExternalRadius; y=-MagnetHeight/2.+MagnetHeightt;};
border b4(t=0., 1.){x=MagnetExternalRadius-(MagnetExternalRadius-MagnetInternalRadius)t; y=MagnetHeight/2.;};
border b5(t=0., 1.){x=MagnetInternalRadius; y=MagnetHeight/2.-MagnetHeightt;};
border b6(t=0., 1.){x=-MagnetInternalRadius-(MagnetExternalRadius-MagnetInternalRadius)t; y=-MagnetHeight/2.;};
border b7(t=0., 1.){x=-MagnetExternalRadius; y=-MagnetHeight/2.+MagnetHeightt;};
border b8(t=0., 1.){x=-MagnetExternalRadius+(MagnetExternalRadius-MagnetInternalRadius)t; y=MagnetHeight/2.;};
border b9(t=0., 1.){x=-MagnetInternalRadius; y=MagnetHeight/2.-MagnetHeightt;};
int nnH = max(2., 10.nnMagnetHeight);
int nnR = max(2., 10.nn(MagnetExternalRadius-MagnetInternalRadius));
mesh Th = buildmesh(
b1(nnL)
+ b2(nnR) + b3(nnH) + b4(nnR) + b5(nnH)
+ b6(-nnR) + b7(-nnH) + b8(-nnR) + b9(-nnH)
);
int Coil1 = Th((MagnetExternalRadius+MagnetInternalRadius)/2., 0.).region;
int Coil2 = Th(-(MagnetExternalRadius+MagnetInternalRadius)/2., 0.).region;
int Box = Th(0., 0.).region;
//Fespace
func Pk = P2;
fespace Ah(Th, Pk);
Ah Az;
//Macro
macro Curl(A) [dy(A), -dx(A)] //
//Current
Ah Jz = J0*(-1.(Th(x,y).region == Coil1) + 1.(Th(x,y).region == Coil2));
plot(Jz, nbiso=30, fill=true, value=true, cmm=“J”, wait=true);
//Problem
func Mu = Mu0 + (MuC-Mu0)*((region==Coil1)+(region==Coil2));
func Nu = 1./Mu;
varf vLaplacian (Az, AAz)
= int2d(Th)(
Nu * Curl(Az)’ * Curl(AAz)
)
- int2d(Th)(
Jz * AAz
)
+ on(BoxWall, Az=0)
;
matrix Laplacian = vLaplacian(Ah, Ah, solver=sparsesolver);
real[int] LaplacianBoundary = vLaplacian(0, Ah);
Az[] = Laplacian^-1 * LaplacianBoundary;
//Magnetic induction
Ah Bx, By;
Bx = dy(Az);
By = -dx(Az);
Ah B = sqrt(Bx^2 + By^2);
//Magnetic field
Ah Hx, Hy;
Hx = Nu * Bx;
Hy = Nu * By;
Ah H = sqrt(Hx^2 + Hy^2);
//Current
Ah J = dx(Hy) - dy(Hx);
//Lorentz force
Ah Lx, Ly;
Lx = -J * By;
Ly = J * Bx;
Ah L = sqrt(Lx^2 + Ly^2);
//Plot
plot(Az, nbiso=20, fill=true, value=true, cmm=“A”, wait=true);
plot(B, nbiso=20, fill=true, value=true, cmm=“B”, wait=true);
plot(Bx, nbiso=20, fill=false, value=true, cmm=“Bx”, wait=true) ;
plot(By, nbiso=20, fill=true, value=true, cmm=“By”, wait=true) ;
plot(H, nbiso=20, fill=true, value=true, cmm=“H”, wait=true);
plot(L, nbiso=20, fill=true, value=true, cmm=“L”, wait=true);
plot(J, nbiso=20, fill=true, value=true, cmm=“Jpost”, wait=true) ;
Hello,
I suspect that there is a wrong sign in the definition of vLaplacian, as
Rot(Rot A)=Rot(B)= mu0.J = Grad(Div A)-Nabla^2(A);
with the (usual) gauge div A=0, the result is Laplacian(A)+mu0*J=0.
If I am right, - int2d(Th)(Jz * AAz) should be +int2d(Th)(Jz * AAz).
HTH
Denis