Boundary conditions for magnetic field simulations

mathematics
1

Hello all,

I’m trying to perform some 2D magnetostatics simulations by solving the poisson equation for the vector potential.

My problem is similar to the magnetostatics example. However in my case I have non negligible currents close to the boundaries, therefore I cannot impose a zero magnetic field boundary condition.

The only option I can come up with to solve this is simulating a much bigger domain. This is not ideal because I have a coupled problem and increasing the size of the domain would make the other parts of the simulations much more expensive.

I have heard of some approach using multiple layers of different materials close to the boundary. Does anyone have any resources on such approach? Any other ideas are also welcome.

Thanks in advance.

1 Like
2

Hello everyone, I was able to solve the problem by myself. The code is not the cleanest and has a lot of things that could be polished but I’ll post it in case anyone encounters a similar problem:

import dolfin as df
import numpy as np
import matplotlib.pyplot as plt

#Mesh
mesh = df.RectangleMesh(df.Point(-10,-10), df.Point(10,10), 50, 50, 'left/right')

#Function Spaces for vector potential and B-Field
V = df.FunctionSpace(mesh, 'CG', 2)
Vdim = df.VectorFunctionSpace(mesh, 'CG', 2)

#Dirichlet BC in all boundaries.
def boundary(x, on_boundary):
      return on_boundary

bc = df.DirichletBC(V, df.Constant(0), boundary)

#Get node coordinates 
bmesh = df.BoundaryMesh(mesh, "exterior", True)

b_coord = bmesh.coordinates()

b_coord_x = np.sort(np.unique(b_coord[:,0]))
b_coord_z = np.sort(np.unique(b_coord[:,1]))

N_layers = 3

#Define meshfunction and set all markers to 0
mshfunc = df.MeshFunction("size_t", mesh, mesh.topology().dim())
mshfunc.set_all(0)

#Overwrite markers of last N_layers cells
for i in reversed(range(N_layers+1)):

    class Layer(df.SubDomain):  
        def inside(self, x, on_boundary):
            return  (df.between(x[0], (b_coord_z[-i-1], b_coord_z[-i])) or df.between(x[1], (b_coord_x[-i-1], b_coord_x[-i]))) or (df.between(x[0], (b_coord_z[i-1], b_coord_z[i])) or df.between(x[1], (b_coord_x[i-1], b_coord_x[i])))
    
    Layer().mark(mshfunc, N_layers+1-i)


# Define integration measure
dx = df.Measure('dx', domain=mesh, subdomain_data=mshfunc)

#Assign permeability to each domain
class Permeability(df.UserExpression): 
    def __init__(self, markers, **kwargs):

        super(Permeability, self).__init__(**kwargs)

        self.markers = markers

    def eval_cell(self, values, x, cell):
        for i in range(N_layers+1):

            if self.markers[cell.index] == i:
                values[0] = (i+1)**2
        
    def value_shape(self):
        return ()

mu = df.interpolate(Permeability(mshfunc), V)

#Currents
class Currents(df.UserExpression):
      def eval(self, value, x):
            if np.sqrt(x[0]**2 + (x[1]-1)**2)< 1:
                  value[0] = 1

            elif np.sqrt(x[0]**2 + (x[1]+1)**2)< 1:
                  value[0] = -1

            else:
                  value[0] = 0


      def value_shape(self):
            return ()

J = df.interpolate(Currents(), V)

#Weak form
A_y = df.TrialFunction(V)

v = df.TestFunction(V)

a =  (1/mu)*df.dot(df.grad(A_y), df.grad(v))*dx

L = J*v*dx

A_y = df.Function(V)

df.solve(a == L, A_y, bc)

#############################
#Postprocessing

x = np.linspace(-10,10, 100)
z = np.linspace(-10,10,100)

X,Z = np.meshgrid(z,x)

psi_grid = [A_y(j,i) for i in x for j in z]
psi_grid = np.reshape(psi_grid, np.shape(X))

mu_grid = [mu(j,i) for i in x for j in z]
mu_grid = np.reshape(mu_grid, np.shape(X))

fig = plt.figure()
ax = fig.add_axes([0,0,1,1])

pcol = ax.pcolor(X, Z, mu_grid)
plt.colorbar(pcol)
cont = ax.contour(X, Z, psi_grid, np.linspace(-0.3, 0.3, 20), colors='cyan', linestyles='solid', linewidths=0.8)

And the outcome is:

bcs
bcs1310×898 186 KB

3

Wait, so what was your solution? You still applied vanishing Dirichlet boundary conditions on the borders of the domain, right? Something you said you could not do, in the 1st post.

4

I implemented a solution similar to what they do in FEMM-Magnetostatics Tutorial.

I am not sure exactly how they implement it there. However, the idea is that when magnetic field traverses an interface between materials with different permeability it deflects (Wikipedia). Then you can use a set of layers of diferent materials in the boundary of the domain and vanishing Dirichlet BCs to simulate an ‘open’ boundary and avoid using a bigger domain than you really need.

Obviously you would need to discard the area where the ‘fictional’ layers are located.

I also have to check how good is the solution really. Anyhow, visually the solution looks good enough. For comparison here are some results solving the same problem for a 5x5 domain with no ‘fictional’ layers, same domain with seven layers, and the same problem solved on a 20x20 domain, respectively.

self_B_0_layers
self_B_7_layers
self_B_20_20

1 Like
5

Interesting approach! I was not aware of the multiple-layer solution.

Did you consider using the coordinate transformation approach? I do not know a place where one could find a detailed description but it is fairly simple and pointed in out in [https://arxiv.org/pdf/2108.05348.pdf] as follows:

“Other problems, such as Poisson’s equation, involve solutions that vary more and more slowly at greater distances—in this case, one can simply employ a coordinate transformation, such as x˜ = tanh x, to remap from x ∈ (−∞, ∞) to x˜ ∈ (−1, 1), and solve the new finite system”

I have tested it in Comsol some time ago with Poisson equation for scalar potential and it worked like a charm. I expect to work for vector potential in cartesian coordinates and most likely in cylindrical coordinates as well. In this way, you only need one fictitious layer except that it is not exactly fictitious but just a compression of the remaining infinity in a finite length of the numerical cell. In principle, one can magnify result of this fictitious space to reconstruct solution in as large a field of view as you like!

2 Likes
6

Thank you Shakeeb, that looks worth considering. I’ll give it a thought and see how it compares to the multimaterial solution.