Hi FEniCS-Community,
I’m very new to FEniCS and I’m trying to simulate the magnetic flux of a cylindrical air coil to get started. I choose this setup because it is possible to compare the results with a Biot-Savart simulation. Later iron will be added to the simulation. The code is similar to the ft11-magnetostatics.py example from the tutorial. I tried to calculate every component of the vector potential independently:
-\nabla \cdot (\mu^{-1} \nabla A_x) = J_x
-\nabla \cdot (\mu^{-1} \nabla A_y) = J_y
-\nabla \cdot (\mu^{-1} \nabla A_z) = J_z
For simplicity I reduced my code here to the calculation of one component. Unfortunately the results are not as expected. The magnetic flux is zero everywhere except the domains where there is a current density. Do I need to add an additional boundary condition? Or is there another problem with my code or my mesh?
I would be very happy if someone finds the mistake/s!
Here you can find the mesh:
The domains were numbered with Gmsh. (Cell index inside the copper of the coil: 4)
I’m using:
Ubuntu 18.04.3
Python 3.6.9
dolfin version 2019.1.0
from fenics import *
from mshr import *
import matplotlib.pyplot as plt
R = 200 #Outer radius air sphere in mm
r = 50 #Radius, coil
J = 1 #Current density in A/m^2
mesh = Mesh("airCoilMesh.xml")
markers = MeshFunction("size_t", mesh,"airCoilMesh_physical_region.xml")
boundaries = MeshFunction("size_t", mesh,"airCoilMesh_facet_region.xml")
# Define function space and integration measure
V = FunctionSpace(mesh, "P", 1)
dx = Measure('dx', domain=mesh, subdomain_data=markers)
dV = Measure("dx")(mesh)
# Define boundary condition
class Boundary(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return x[0]**2+x[1]**2+x[2]**2 >= R**2-tol
boundary = Boundary()
bc = DirichletBC(V, Constant(0), boundary)
# define current densities as function of subdomains
class CurrentDensity_x(UserExpression):
def __init__(self, markers, **kwargs):
super().__init__(**kwargs)
self.markers = markers
def eval_cell(self, values, x, cell):
if (self.markers[cell.index] == 4):
values[0] = -x[1]/r*J
else:
values[0] = Constant(0.0)
# permeability
class Permeability(UserExpression):
def __init__(self, markers, **kwargs):
super().__init__(**kwargs)
self.markers = markers
def eval_cell(self, values, x, cell):
if self.markers[cell.index] == 4:
values[0] = -6.4e-6 # copper
else:
values[0] = 4*pi*1e-7 # vacuum
mu = Permeability(markers, degree=1)
# initialise current densities
JX = CurrentDensity_x(markers, degree = 1)
Jx = project(JX, V)
# define variational problems for x, y, z
Ax = TrialFunction(V)
vx = TestFunction(V)
# left hand sides of variational problems
ax = dot(grad(Ax), grad(vx))*dV
# right hand sides of variational problems
Lx = mu*Jx*vx*dV
# solve variational problems for magnetic Ax (in Vs/mm)
Ax = Function(V)
solve(ax == Lx, Ax, bc)
# plot solutions
plt.figure()
plot(Ax)
plt.show()