Magnetic field in 3D of a "permanent magnet"

[raboynics]

I am trying to calculate the magnetic field of a permanent magnet in 3D (just for fun), in the region outside of the magnet (because I have read that in general a permanent magnet is ferromagnetic and the relation between its magnetization and the field is non linear, i.e. very complicated, which may matter if more than 1 magnet is considered at once).

I took the magnet as spherical and the surrounding space as being a box. Here is my gmsh file:

SetFactory(“OpenCASCADE”);
Mesh.RandomFactor=1.0e-4;
Sphere(1) = {0, 0, 0, 1};
Box(2) = {-4.9, -4.9, -4.9, 10, 10, 10};
BooleanDifference{ Volume{2}; Delete; }{ Volume{1}; Delete; }
Physical Volume(“testing_vol”) = {2};
Physical Surface(“sphere_surface_magnet”) = {1};
Physical Surface(“box_infinity_surfaces”) = {2, 6, 5, 4, 7, 3};

I generated the mesh with gmsh GUI and refined it there too. So it looks like the following picture

Even though it isn’t visible, there’s really a spherical empty space near the middle of the box, the region of the permanent magnet and where I am not interested to calculate the magnetic field.

I have then followed the magnetostatics FEniCS tutorial (https://fenicsproject.org/pub/tutorial/html/._ftut1015.html) and tried to extend it to 3D. In my case the \vec B field isn’t two dimensional, it is 3D. And the vector potential isn’t along z only (it could have a spherical symmetry because I picked a spherical magnet but I would rather not enforce such symmetry since I plan to modify the geometry. Furthermore it isn’t exactly centered in the box, so the solution shouldn’t be exactly symmetric with respect to rotations).
My FEniCS code (which works, i.e. returns something) is the following

 import numpy as np
 import meshio
 msh = meshio.read("sphere_magnet_test.msh")
 
 ''' # Uncomment to create the files for the mesh to use with FEniCS.
 cells = np.vstack(np.array([cells.data for cells in msh.cells if cells.type == "tetra"]))
 tetra_data = msh.cell_data_dict["gmsh:physical"]["tetra"]
 tetra_mesh = meshio.Mesh(points=msh.points, cells=[("tetra", cells)], cell_data={"name_to_read": [tetra_data]})
 facet_cells = np.vstack([cells.data for cells in msh.cells if cells.type == "triangle"])
 facet_data = msh.cell_data_dict["gmsh:physical"]["triangle"]
 facet_mesh = meshio.Mesh(points=msh.points,
                          cells=[("triangle", facet_cells)],
                          cell_data={"name_to_read": [facet_data]})
 # Write mesh
 meshio.xdmf.write("mesh_test0.xdmf", tetra_mesh)
 meshio.xdmf.write("facet_mesh_test0.xdmf", facet_mesh)
 '''
 #exit()
 
 from dolfin import *
 mesh = Mesh()
 
 # Read mesh.
 with XDMFFile("mesh_test0.xdmf") as infile:
     infile.read(mesh)
 mvc = MeshValueCollection("size_t", mesh, 2)
 with XDMFFile("facet_mesh_test0.xdmf") as infile:
     infile.read(mvc, "name_to_read")
 mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
 mvc2 = MeshValueCollection("size_t", mesh, 3)
 with XDMFFile("mesh_test0.xdmf") as infile:
     infile.read(mvc2, "name_to_read")
 cf = cpp.mesh.MeshFunctionSizet(mesh, mvc2)
 dx = Measure("dx", domain=mesh,subdomain_data=cf)
 ds = Measure("ds", domain=mesh, subdomain_data=mf)
 dS = Measure("dS", domain=mesh, subdomain_data=mf)
 subdomains = cf
 boundary_parts = mf
 
 print('Expected volume:', 10**3-(4.0/3)*np.pi*1**3)
 print('Computed volume: ',assemble(Constant(1)*dx(domain=mesh, subdomain_data=cf)))
 print('computed volume again: ',assemble(Constant(1)*dx(domain=mesh, subdomain_data=cf, subdomain_id=1)))
 
 print('Expected magnet surface:', 4*np.pi*1**2)
 print('Magnet surface: ',assemble(Constant(1)*ds(domain=mesh, subdomain_data=mf, subdomain_id=2)))
 
 print('Box expected surface: ', 10**2*6)
 print('Box computed surface: ', assemble(Constant(1)*ds(domain=mesh, subdomain_data=mf, subdomain_id=3)))
 
 # Define function space
 V = FunctionSpace(mesh, 'P', 1)
 
 # Define parameters
 #magnetic_potential_sphere = Constant(0)
 magnetic_potential_walls = Constant(0)
 
 effective_current = Constant(1.0)
 
 # Vacuum permeability.
 mu_0 = 4*np.pi*1e-7
 
 #bc0 = DirichletBC(V, magnetic_potential_walls, mf, 2)
 bc1 = DirichletBC(V, magnetic_potential_walls, mf, 3)
 #bcs = [bc0, bc1]
 
 # Define variational problem
 A_z = TrialFunction(V) # This doesn't work when the eq. is non linear.
 A_x = TrialFunction(V)
 A_y = TrialFunction(V)
 
 
 v = TestFunction(V)
 a = (1 / mu_0)*dot(grad(A_z), grad(v))*dx
 L = effective_current * v * dx
 
 a1 = (1 / mu_0)*dot(grad(A_x), grad(v))*dx
 L1 = effective_current * v * dx
 
 a2 = (1 / mu_0)*dot(grad(A_y), grad(v))*dx
 L2 = effective_current * v * dx
 
 
 # Solve variational problem
 A_z = Function(V)
 solve(a == L, A_z, bcs = bc1)
 
 A_x = Function(V)
 solve(a1 == L1, A_x, bcs = bc1)
 
 A_y = Function(V)
 solve(a2 == L2, A_y, bcs = bc1)
 
 # Compute magnetic field (B = curl A)
 W = VectorFunctionSpace(mesh, 'P', 1)
 B = project(as_vector( (A_z.dx(1)-A_y.dx(2), -A_z.dx(0)+A_x.dx(2), A_y.dx(0)-A_x.dx(1)) ), W)
 
 # Save solution to file
 vtkfile_A_z = File('magnetostatics_test/potential_z.pvd')
 vtkfile_A_x = File('magnetostatics_test/potential_x.pvd')
 vtkfile_A_y = File('magnetostatics_test/potential_y.pvd')
 vtkfile_B = File('magnetostatics_test/field.pvd')
 vtkfile_A_z << A_z
 vtkfile_A_x << A_x
 vtkfile_A_y << A_y
 vtkfile_B << B

As you can see, since the potential vector is really a 3D vector in my case, I solve 3 variational problems, one for each component of \vec A. Is it a correct approach? Or should have I used a mixed element approach instead? Or simply just a single variational problem without ME?

I would have expected the \vec B field to look like a dipole (just like the Earth’s magnetic field), but it doesn’t look like that, to me. Here’s a view of the field with Paraview:

To me, it looks like \vec B is circular around an axis. I probably picked wrong boundary conditions for the vector potential. If I assume that the magnet is polarized in a particular direction then I expect the boundary condition on \vec A to be different according the direction, which is not my case. In my case it’s as if the polarization of the magnet had the same x, y and z components, maybe explaining why the \vec B field “circulates” around an axis. I am not quite sure where I go wrong. I also assumed that the magnet’s polarization is equivalent to an effective current, i.e. as if the \vec B field was produced by a current density rather than by a magnetic material.

Edit: By modyfing a little bit the code, more specifically I changed a part for this one instead:

magnetic_potential_walls = Constant(0)
magnetic_potential_sphere = Constant(1000000)

effective_current_z = Constant(0.0)
#effective_current_z = Expression('10000*atan(sqrt(x[0]*x[0]+x[1]*x[1])/x[2])',degree=2)

effective_current_x = Constant(0.0)
effective_current_y = Constant(0.0)

# Magnetic permeability.
mu_0 = 4*np.pi*1e-7

bc_walls = DirichletBC(V, magnetic_potential_walls, mf, 3)
bc_sphere_z = DirichletBC(V, magnetic_potential_sphere, mf, 2)
bc_sphere_x = DirichletBC(V, magnetic_potential_walls, mf, 2)
bc_sphere_y = DirichletBC(V, magnetic_potential_walls, mf, 2)

# Define variational problem
A_z = TrialFunction(V) # This doesn't work when the eq. is non linear.
#A_z = Function(V) # works for non linear equations.
A_x = TrialFunction(V)
A_y = TrialFunction(V)


v = TestFunction(V)

a = (1 / mu_0)*dot(grad(A_z), grad(v))*dx
L = effective_current_z * v * dx

a1 = (1 / mu_0)*dot(grad(A_x), grad(v))*dx
L1 = effective_current_x * v * dx

a2 = (1 / mu_0)*dot(grad(A_y), grad(v))*dx
L2 = effective_current_y * v * dx


# Solve variational problem
A_z = Function(V)
solve(a == L, A_z, bcs = [bc_sphere_z, bc_walls])

A_x = Function(V)
solve(a1 == L1, A_x, bcs = [bc_sphere_x, bc_walls])

A_y = Function(V)
solve(a2 == L2, A_y, bcs = [bc_sphere_y, bc_walls])

So that there’s no electric current at all anymore, only different values of the potential vector on the sphere vs on the box walls. I specified the anisotropy in A_z via the potential vector boundary condition on the sphere’s surface.
I get the following magnetic field:

Unfortunately it doesn’t look correct to me. The \vec B field should point either towards or outwards the sphere in the region closest to the sphere, not around it.

[kamensky]

The linked example’s formulation for A_z is based on the assumption that fields are constant in the z-direction, and that the magnetic field is entirely in-plane. This is not the case in your example, and it is not valid to make this assumption for each component of \mathbf{A}, then superpose the results.

The standard approach for solving 3D magnetostatic problems is to solve the “curl-curl” problem for the vector potential \mathbf{A}, discretizing it with a special type of vector-valued finite element, called a Nédélec element. Nédélec function spaces are subsets of the space H(\text{curl}) of functions whose curls are square-integrable. However, unlike standard C^0 finite element spaces, Nédélec spaces are not subsets of H^1 (i.e., their full gradients are not square-integrable, only their curls). The practical implication of this is that tangential components of Nédélec functions are continuous across element faces, while normal components are not.

Clearly H^1-conforming functions are a subset of H(\text{curl}), but are not in general suitable for the curl-curl problem, because this subset is too restricted, as illustrated using FEniCS here, in the context of eigenvalue problems:

https://fenics.readthedocs.io/projects/dolfin/en/2017.2.0/demos/maxwell-eigenvalues/python/demo_maxwell-eigenvalues.py.html

For further discussion on subtleties of setting up curl-curl problems in FEniCS, you can take a look at this old discussion:

https://fenicsproject.org/qa/11451/nedelec-conforming-elements-converge-simple-test-problem/

[raboynics]

Thank you very much for replying!
This is above my head at the moment, I will try to see how to proceed.

However before I dive into it, I do not understand why the problem isn’t decoupled as I thought it is.
From Maxwell’s equation \nabla \times \vec H = \frac{\partial \vec D}{\partial t}+\vec J with \vec D=\vec 0 and assuming a linear relation between \vec H and \vec B (valid in vacuum or air, which is what matters to me, since this is where I solve the problem), one finds that \nabla \times \left( \frac{1}{\mu} \nabla \times \vec A \right)=\vec J . If \nabla \mu = \vec 0 (usually assumed to hold), then we get \nabla (\nabla \cdot \vec A)-\nabla ^2 \vec A=\mu \vec J where \nabla^2 is the “vector Laplacian” operator (this is eq.5.92 in Jackson’s Classical Electrodynamics textbook), so that \nabla ^2 \vec A= (\nabla^2 A_x, \nabla^2 A_y, \nabla^2 A_z). So this is a vector equation, to me it looks like every component is decoupled to any other. So solving this equation, to me, looks like to solve the equation for each component.

To me this leads to solve (assuming the Coulomb gauge \nabla \cdot \vec A =0)
\nabla^2 A_x=-\mu J_x for A_x (this solves the first component of the vector equation),

\nabla^2 A_y=-\mu J_y for A_y (this solves the second component of the vector equation),

and \nabla^2 A_z=-\mu J_z for A_z (this solves the third and last component of the vector equation).

Here, in these 3 equations, \nabla^2 is the (scalar) Laplacian.

I do not see where I go wrong. Could you please comment?

[kamensky]

If I understand your argument correctly, the logic is that \nabla\cdot\mathbf{A} = 0 implies that each component of \mathbf{A} satisfies an independent Poisson equation. However, if you solve three independent Poisson equations, I don’t see why the resulting vector field would be guaranteed to be solenoidal. Thus, the components would still be coupled through the constraint that \nabla\cdot\mathbf{A} = 0.

[nate]

For details consider Peter Monk’s excellent book, or from an engineering perspective, Jin’s excellent book.

[raboynics]

Thanks for the references, but I am not sure what I should look for, exactly. For instance in Jin’s textbook, it is mentioned that for 3D magnetostatics problems, it is indeed question of solving for the vector potential (and not the “curl curl” way mentioned by Kamensky). He mentions which boundary conditions to implement, but I have yet to get familiar with his notation. I suspect that if I could fix/tweak the boundary conditions in my code, it should work (i.e. return the correct result). But overall, I am still missing what I do wrong, exactly.

Regarding the first reference, I do not see anything related to magnetostatics in particular in the book. It has a part dealing with vector potentials in general, but I fail to see the relation with Maxwell’s equations. The “Coulomb” gauge for example is not mentioned (in neither book for that matter, though it is at least described in Jin’s).

So, yeah, I still do not see what is wrong with my approach. I suspect it’s just the boundary conditions, but that’s a guess.

[raboynics]

Essentially yes. This is written explicitely in Jackson’s textbook (p.181 of the 3rd edition). It should be Poisson equation (since the current density is not necessarily null), but in the case of no electric current this is indeed a Laplace equation for each of the components of the vector potential. I am not sure yet how to answer to your remark regarding the solenoidal condition. Good point to ponder on.