Magnetostatic issue

dega · November 4, 2022, 11:16am

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 copperRadius   = 8e-3;
real CoilYSpacing   = 0.018;
real CoilRadius     = 0.015;

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

int BoxWall = 10;

int nnL = max(2., 5*BoxRadius*nn/2.);
border b1(t=0., 2.*pi){x=BoxRadius*cos(t); y=BoxRadius*sin(t); label=BoxWall;};

border b2(t=0., 2.*pi){x=CoilRadius+copperRadius*cos(t);  y=-(CoilYSpacing*2.)+copperRadius*sin(t);};
border b3(t=0., 2.*pi){x=CoilRadius+copperRadius*cos(t);  y=-CoilYSpacing+copperRadius*sin(t);};
border b4(t=0., 2.*pi){x=CoilRadius+copperRadius*cos(t);  y=copperRadius*sin(t);};
border b5(t=0., 2.*pi){x=CoilRadius+copperRadius*cos(t);  y=CoilYSpacing+copperRadius*sin(t);};
border b6(t=0., 2.*pi){x=CoilRadius+copperRadius*cos(t);  y=(CoilYSpacing*2.)+copperRadius*sin(t);};

border b7(t=0., 2.*pi){x=-CoilRadius+copperRadius*cos(t); y=(-CoilYSpacing*2.)+copperRadius*sin(t);};
border b8(t=0., 2.*pi){x=-CoilRadius+copperRadius*cos(t); y=-CoilYSpacing+copperRadius*sin(t);};
border b9(t=0., 2.*pi){x=-CoilRadius+copperRadius*cos(t); y=copperRadius*sin(t);};
border b10(t=0., 2.*pi){x=-CoilRadius+copperRadius*cos(t); y=CoilYSpacing+copperRadius*sin(t);};
border b11(t=0., 2.*pi){x=-CoilRadius+copperRadius*cos(t); y=(CoilYSpacing*2.)+copperRadius*sin(t);};


int nnC = max(2., 900*copperRadius*nn/2.);

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 Coil1 = Th(CoilRadius, -2.*CoilYSpacing).region;
int Coil2 = Th(CoilRadius, -CoilYSpacing).region;
int Coil3 = Th(CoilRadius, 0.).region;
int Coil4 = Th(CoilRadius, CoilYSpacing).region;
int Coil5 = Th(CoilRadius, 2.*CoilYSpacing).region;

int Coil6 = Th(-CoilRadius, -2.*CoilYSpacing).region;
int Coil7 = Th(-CoilRadius, -CoilYSpacing).region;
int Coil8 = Th(-CoilRadius, 0.).region;
int Coil9 = Th(-CoilRadius, CoilYSpacing).region;
int Coil10 = Th(-CoilRadius, 2.*CoilYSpacing).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) - 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");

dega · November 4, 2022, 11:19am

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 copperRadius   = 8e-3;
real CoilYSpacing   = 0.018;
real CoilRadius     = 0.015;

int BoxWall = 10;

int nnL = max(2., 5*BoxRadius*nn/2.);
border b1(t=0., 2.*pi){x=BoxRadius*cos(t); y=BoxRadius*sin(t); label=BoxWall;};

border b2(t=0., 2.*pi){x=CoilRadius+copperRadius*cos(t);  y=-CoilYSpacing+copperRadius*sin(t);};
border b3(t=0., 2.*pi){x=CoilRadius+copperRadius*cos(t);  y=copperRadius*sin(t);};
border b4(t=0., 2.*pi){x=CoilRadius+copperRadius*cos(t);  y=CoilYSpacing+copperRadius*sin(t);};
border b5(t=0., 2.*pi){x=CoilRadius+copperRadius*cos(t);  y=(CoilYSpacing*2.)+copperRadius*sin(t);};

border b6(t=0., 2.*pi){x=-CoilRadius+copperRadius*cos(t); y=-CoilYSpacing+copperRadius*sin(t);};
border b7(t=0., 2.*pi){x=-CoilRadius+copperRadius*cos(t); y=copperRadius*sin(t);};
border b8(t=0., 2.*pi){x=-CoilRadius+copperRadius*cos(t); y=CoilYSpacing+copperRadius*sin(t);};
border b9(t=0., 2.*pi){x=-CoilRadius+copperRadius*cos(t); y=(CoilYSpacing*2.)+copperRadius*sin(t);};


int nnC = max(2., 900*copperRadius*nn/2.);

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 Coil1 = Th(CoilRadius, -CoilYSpacing).region;
int Coil2 = Th(CoilRadius, 0.).region;
int Coil3 = Th(CoilRadius, CoilYSpacing).region;
int Coil4 = Th(CoilRadius, 2.*CoilYSpacing).region;

int Coil5 = Th(-CoilRadius, -CoilYSpacing).region;
int Coil6 = Th(-CoilRadius, 0.).region;
int Coil7 = Th(-CoilRadius, CoilYSpacing).region;
int Coil8 = Th(-CoilRadius, 2.*CoilYSpacing).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) - 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");

dega · November 4, 2022, 1:11pm

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

dega · November 4, 2022, 1:12pm

The 4 coil one seems to have computed ok:

gero · November 5, 2022, 9:36pm

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

dega · November 6, 2022, 10:33am

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

gero · November 8, 2022, 12:44pm

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

dega · November 8, 2022, 1:31pm

Ah I see…
Defnitely my mistake and no bug then!

Thanks again,
Marco

Topic		Replies	Views
Cyclic boundary conditions General Discussion	0	288	March 26, 2021
Improving residual error, magnetostatic General Discussion	2	226	July 7, 2022
Magnetostatic Simulation of Dipole General Discussion	1	879	March 7, 2022
Magnetostatics inconsistent sign for current and field General Discussion	3	629	January 14, 2020
Lines missing in 3D algorithm of magnetostatic module General Discussion	7	480	May 3, 2022

Magnetostatic issue

Related Topics