### Vector potential and magnetic field calculated as zero

97
views
0
4 weeks ago by

Hello everyone,

I'm trying to solve a linear magnetic field problem. My code so far is this:

# some constants
mu_ri = 1.0e6
mu_rc = 1.0
n = 50

# coil currents
I_left = 10.0
I_right = 10.0

# following current densities
J_left = I_left*n/(109.4*(27.88))
J_right = I_right*n/(109.4*(27.88))

mesh = Mesh("Feldgenerator.xml")
subdomains = MeshFunction("size_t", mesh, "Feldgenerator_physical_region.xml")

# define functions
V = VectorFunctionSpace(mesh, "CG", 2)
A = TrialFunction(V)
v = TestFunction(V)

# define boundary conditions
bcs = DirichletBC(V, (0.0, 0.0, 0.0), 'on_boundary')
dx = Measure("dx")(subdomain_data = subdomains)

# define relative permeability for subdomains
class Permeability(Expression):
def __init__(self, subdomains, **kwargs):
self.subdomains = subdomains
def eval_cell(self, values, x, cell):
if (self.subdomains[cell.index] == 0) or (self.subdomains[cell.index] == 18) or (self.subdomains[cell.index] == 19):
values[0] = 1.0                                 # 0 background, 18/19 isolations
elif self.subdomains[cell.index] == 1:
values[0] = mu_ri                              # iron yoke
else:
values[0] = mu_rc

# permeability
mu = 4*pi*1e-7*Permeability(subdomains, degree = 1)

# left hand side of variational problem

# define current density vectorsand  right hand side of variational problem
# left coil
Jl1 = project(as_vector((J_left, 0.0, 0.0)), V)                             # top
Ll1 = mu*inner(Jl1, v)*dx(2)
Jl2 = project(as_vector((1/sqrt(2)*J_left, 0.0, -1/sqrt(2)*J_left)), V)     # top left
Ll2 = mu*inner(Jl2, v)*dx(3)
Jl3 = project(as_vector((0.0, 0.0, -J_left)), V)                            # left
Ll3 = mu*inner(Jl3, v)*dx(4)
Jl4 = project(as_vector((-1/sqrt(2)*J_left, 0.0, -1/sqrt(2)*J_left)), V)    # bottom left
Ll4 = mu*inner(Jl4, v)*dx(5)
Jl5 = project(as_vector((-J_left, 0.0, 0.0)), V)                            # bottom
Ll5 = mu*inner(Jl5, v)*dx(6)
Jl6 = project(as_vector((-1/sqrt(2)*J_left, 0.0, 1/sqrt(2)*J_left)), V)     # bottom right
Ll6 = mu*inner(Jl6, v)*dx(7)
Jl7 = project(as_vector((0.0, 0.0, J_left)), V)                             # right
Ll7 = mu*inner(Jl7, v)*dx(8)
Jl8 = project(as_vector((1/sqrt(2)*J_left, 0.0, 1/sqrt(2)*J_left)), V)      # top right
Ll8 = mu*inner(Jl8, v)*dx(9)

Ll = Ll1 + Ll2 + Ll3 + Ll4 + Ll5 + Ll6 + Ll7 + Ll8

# right coil
Jr1 = project(as_vector((-J_right, 0.0, 0.0)), V)                            # top
Lr1 = mu*inner(Jl1, v)*dx(10)
Jr2 = project(as_vector((-1/sqrt(2)*J_right, 0.0, -1/sqrt(2)*J_right)), V)   # top left
Lr2 = mu*inner(Jl2, v)*dx(11)
Jr3 = project(as_vector((0.0, 0.0, -J_right)), V)                            # left
Lr3 = mu*inner(Jl3, v)*dx(12)
Jr4 = project(as_vector((1/sqrt(2)*J_right, 0.0, -1/sqrt(2)*J_right)), V)    # bottom left
Lr4 = mu*inner(Jl4, v)*dx(13)
Jr5 = project(as_vector((J_right, 0.0, 0.0)), V)                             # bottom
Lr5 = mu*inner(Jl5, v)*dx(14)
Jr6 = project(as_vector((1/sqrt(2)*J_right, 0.0, 1/sqrt(2)*J_right)), V)     # bottom right
Lr6 = mu*inner(Jl6, v)*dx(15)
Jr7 = project(as_vector((0.0, 0.0, J_right)), V)                             # right
Lr7 = mu*inner(Jl7, v)*dx(16)
Jr8 = project(as_vector((-1/sqrt(2)*J_right, 0.0, 1/sqrt(2)*J_right)), V)    # top right
Lr8 = mu*inner(Jl8, v)*dx(17)

Lr = Lr1 + Lr2 + Lr3 + Lr4 + Lr5 + Lr6 + Lr7 + Lr8

L = Ll + Lr

aa = assemble(a)
LL = assemble(L)
bcs.apply(aa,LL)

# solve variational problem
A = Function(V)                 # in Vs/mm
solve(aa, A.vector(), LL)

Ax, Ay, Az = A.split(deepcopy=True)

# compute magnetic field (B = curl(A)), in Vs/m^2 = T
Bx = (Az.dx(1) - Ay.dx(2))*1.0e6
By = (Ax.dx(2) - Az.dx(0))*1.0e6
Bz = (Ay.dx(0) - Ax.dx(1))*1.0e6
B = project(as_vector((Bx, By, Bz)), V)

vtkfile_A = File("solutions/potential.pvd")
vtkfile_A << A

vtkfile_B =  File("solutions/field.pvd")
vtkfile_B << B

# print values of magnetic field components on vertices
for v in vertices(mesh):
x = v.point().x()
y = v.point().y()
z = v.point().z()
print(x, y, z, Bx(x, y, z), By(x, y, z), Bz(x, y, z))

# print magnetic field
plot(B, mode = "glyphs")
plt.show()

This is a geometry consisting of an iron yoke surrounded by air, that has two coils at its ends. These coils are separated into various subdomains in order to define the current densities J more accurately, that's why the calculation of the J's and L's is so long.

Right now, I'm stuck with two problems:

1. When I try to run this code, I get an error message with the print function. It says "AttributeError: 'float' object has no attribute 'get'".What does that mean? I guess it is because I'm not defining the vector/function B and its components properly, but I'm not sure how to fix it. I had it working half an hour ago - just couldn't identify what changed... I tried omitting the factor *1.0e6 when calculating the B components, same error still.

2. The more general problem is with the solution I get without the printing error. The values for both the vector potential A and the magnetic field B are zero everywhere, except inside the coils. Considering the physics this doesn't make any sense at all. I was unsure if that is because the right-hand side of the variational problem, L, is non-zero only in the coils (because J=0 everywhere else), but that should be correct, shouldn't it? Is there a mistake in the definition of my variational problem?

Btw, for both problems I get the same result when I just compute B as

B = curl(A)

I'm totally new to FEniCS, so I'd really appreciate any suggestions for correcting this stuff.

Rebecca

Community: FEniCS Project

1
4 weeks ago by
Hi,

your problem 1 is (probably) because of
Bx(x, y, z)
you need to "calculate" the nodal values first by using a project, and then you can use such a feature. For example, have a look at your magnetic flux,
B = project(as_vector((Bx, By, Bz)), V)
this is now a function defined everywhere (the nodal values are stored, in between the form functions defined in V are used for computing the values of B vector). Now you can ask the value in the x,y,z loop like that
B( Point(x,y,z) )[0] for x-component
B( Point(x,y,z) )[1] for y-component
B( Point(x,y,z) )[2] for z-component

Your problem 2 is (most likely) because of the boundary conditions. If you set far away conditions as 0 then the solution is zero as well (unless the body brings in a disturbance, within the coil). From the physical point of view, your formulation does not make sense, I believe that you want to solve an alternative current, otherwise, you will not induce a magnetic field. But this forum is for FEniCS specific questions not for physics related.

Best, Emek

Hi Emek,

1. Perfect, that's exactly what I was looking for. Working fine now!

2. Unfortunately, the physics are correct in this case. A constant current induces a constant magnetic field, which is what I need. The PDE I'm solving is Laplace(A) = -mu*J with a time-constant J. In a simple 2D simulation, this worked perfectly fine, also with the zero boundary conditions. So it has to be something with the code... I took your advice and had a look at the bcs. I tried omitting them, same result still. Also when I set them to 100, the result was the same - zero even on the boundaries. That shouldn't be, should it?

Anyway, thanks again for your help.

written 4 weeks ago by Rebecca Pflanze
This equation comes from the Maxwell equation, after using Lorenz gauge and neglecting second time derivative of A, known as steady state. If you assume steady state, then you expect from the Newton solver to find a solution far far away from the zero solution (initial conditions). Your physical assumptions are very risky for numerical reasons. They are valid but only in very very special purposes, in textbooks you find this assumption in order to construct a problem with a closed form solution. In numerical simulation, you do not need these assumptions.

Once more, if you want to simulate coils, transformers and so on, you need AC not DC, otherwise, you will not induce a magnetic field (generating the magnetic flux).
written 4 weeks ago by Emek