Magnetostatic issue

Hello,
I am trying to make a simple magnetostatic simulation of a coil, I started from the example in:
https://modules.freefem.org/modules/magnetostatic/

I created a simple gemometry with 4 coils and the results seem ok, when I add one more, the simulation seems to break

This is the code with the problem:

``````//Parameters
real Mu0 = 4.*pi*1.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 CoilYSpacing   = 0.018;

real pipeHeight     =0.3;
real pipeThickness  =0.01;
real pipeID         =0.04;

int BoxWall = 10;

mesh Th = buildmesh(
b1(nnL) + b2(nnC) + b3(nnC) + b4(nnC) + b5(nnC) + b6(nnC) + b7(nnC) + b8(nnC) + b9(nnC) + b10(nnC) + b11(nnC)
);

plot(Th,wait=1);

//cout << b13;

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) - 1.*(Th(x,y).region == Coil3) - 1.*(Th(x,y).region == Coil4) - 1.*(Th(x,y).region == Coil5)
+1.*(Th(x,y).region == Coil6) + 1.*(Th(x,y).region == Coil7) + 1.*(Th(x,y).region == Coil8) + 1.*(Th(x,y).region == Coil9) + 1.*(Th(x,y).region == Coil10));
plot(Jz, nbiso=30, fill=true, value=true, cmm="J - Current Density");

//Problem
func Mu = Mu0 + (MuC-Mu0)*((region==Coil1)+(region==Coil2)+(region==Coil3)+(region==Coil4)+(region==Coil5)+(region==Coil6)+(region==Coil7)+(region==Coil8)+(region==Coil9)+(region==Coil10));

func Nu = 1./Mu;
varf vLaplacian (Az, AAz)
= int2d(Th)(
Nu * Curl(Az)' * Curl(AAz)
)
- int2d(Th)(
Jz * AAz
)
+ on(BoxWall, Az=0)
;
matrix<real> 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=30, fill=true, value=true, cmm="A");
plot(B, nbiso=30, fill=true, value=true, cmm="B - Magnetic induction");
plot(H, nbiso=30, fill=true, value=true, cmm="H - Magnetic field");
plot(J, nbiso=30, fill=true, value=true, cmm="J - Current");
plot(L, nbiso=30, fill=true, value=true, cmm="Lorentz force");
Ah MuVal =  Mu;
plot(MuVal, nbiso=30, fill=true, value=true, cmm="Permeability");

``````
1 Like

For reference,
this is the code with just 4 coils, that doesnâ€™t seem to have the issue:

``````//Parameters
real Mu0 = 4.*pi*1.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 CoilYSpacing   = 0.018;

int BoxWall = 10;

mesh Th = buildmesh(
b1(nnL) + b2(nnC) + b3(nnC) + b4(nnC) + b5(nnC) + b6(nnC) + b7(nnC) + b8(nnC) + b9(nnC)
);

plot(Th,wait=1);

//cout << b13;

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) - 1.*(Th(x,y).region == Coil3) - 1.*(Th(x,y).region == Coil4)
+1.*(Th(x,y).region == Coil5) + 1.*(Th(x,y).region == Coil6) + 1.*(Th(x,y).region == Coil7) + 1.*(Th(x,y).region == Coil8));
plot(Jz, nbiso=30, fill=true, value=true, cmm="J - Current Density");

//Problem
func Mu = Mu0 + (MuC-Mu0)*((region==Coil1)+(region==Coil2)+(region==Coil3)+(region==Coil4)+(region==Coil5)+(region==Coil6)+(region==Coil7)+(region==Coil8));

func Nu = 1./Mu;
varf vLaplacian (Az, AAz)
= int2d(Th)(
Nu * Curl(Az)' * Curl(AAz)
)
- int2d(Th)(
Jz * AAz
)
+ on(BoxWall, Az=0)
;
matrix<real> 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=30, fill=true, value=true, cmm="A");
plot(B, nbiso=30, fill=true, value=true, cmm="B - Magnetic induction");
plot(H, nbiso=30, fill=true, value=true, cmm="H - Magnetic field");
plot(J, nbiso=30, fill=true, value=true, cmm="J - Current");
plot(L, nbiso=30, fill=true, value=true, cmm="Lorentz force");
Ah MuVal =  Mu;
plot(MuVal, nbiso=30, fill=true, value=true, cmm="Permeability");

``````

A screenshot to show the odd solution for the 5 coil:

The 4 coil one seems to have computed ok:

I donâ€™t know if this helps:

I could not figure out what is wrong with your codeâ€¦ so I reworked the script to automatically generate the mesh for multiple windings (via nWinding; using array of border {see documentation}) and to reassign region numbers (as Frederic Hecht suggested in a recent answer to a different issue).

It seems to work for nWinding 3â€¦7. At least the results look symmetric.

revN5.edp (2.9 KB)

Perhaps it helps you or somebody else to figure out what goes wrong with your scripts.

Kind regards
Gerald

1 Like

Hi Gero,
I like your code, way neater than my raw one.

Maybe it has somethign to do with region names, a small bug lurking somewhereâ€¦

Thank you very much!
Marco

1 Like

Just to conclude the matter:

The original problem is caused by a poor choice of the value of BoxWall
int BoxWall = 10;
as for higher numbers of windings 10 is also used for an automatically assigned label.

Thus the boundary condition is applied incorrectly.

Lesson learned:
Use reasonably high numbers for manual assigned labelsâ€¦
so e.g.
int BoxWall = 99;
fixes the problemâ€¦

All the best
Gerald

Ah I seeâ€¦
Defnitely my mistake and no bug then!

Thanks again,
Marco