Permalink
Find file
a516f1b Nov 24, 2016
104 lines (81 sloc) 2.82 KB
"""
FEniCS tutorial demo program: Magnetic field generated by a copper
wire wound around an iron cylinder. The solution is computed by
solving the Poisson equation for the z-component of the magnetic
vector potential.
-Laplace(A_z) = mu_0 * J_z
"""
from __future__ import print_function
from fenics import *
from mshr import *
from math import sin, cos, pi
a = 1.0 # inner radius of iron cylinder
b = 1.2 # outer radius of iron cylinder
c_1 = 0.8 # radius for inner circle of copper wires
c_2 = 1.4 # radius for outer circle of copper wires
r = 0.1 # radius of copper wires
R = 5.0 # radius of domain
n = 10 # number of windings
# Define geometry for background
domain = Circle(Point(0, 0), R)
# Define geometry for iron cylinder
cylinder = Circle(Point(0, 0), b) - Circle(Point(0, 0), a)
# Define geometry for wires (N = North (up), S = South (down))
angles_N = [i*2*pi/n for i in range(n)]
angles_S = [(i + 0.5)*2*pi/n for i in range(n)]
wires_N = [Circle(Point(c_1*cos(v), c_1*sin(v)), r) for v in angles_N]
wires_S = [Circle(Point(c_2*cos(v), c_2*sin(v)), r) for v in angles_S]
# Set subdomain for iron cylinder
domain.set_subdomain(1, cylinder)
# Set subdomains for wires
for (i, wire) in enumerate(wires_N):
domain.set_subdomain(2 + i, wire)
for (i, wire) in enumerate(wires_S):
domain.set_subdomain(2 + n + i, wire)
# Create mesh
mesh = generate_mesh(domain, 32)
# Define function space
V = FunctionSpace(mesh, 'P', 1)
# Define boundary condition
bc = DirichletBC(V, Constant(0), 'on_boundary')
# Define subdomain markers and integration measure
markers = MeshFunction('size_t', mesh, 2, mesh.domains())
dx = Measure('dx', domain=mesh, subdomain_data=markers)
# Define current densities
J_N = Constant(1.0)
J_S = Constant(-1.0)
# Define magnetic permeability
class Permeability(Expression):
def __init__(self, mesh, **kwargs):
self.markers = markers
def eval_cell(self, values, x, cell):
if markers[cell.index] == 0:
values[0] = 4*pi*1e-7 # vacuum
elif markers[cell.index] == 1:
values[0] = 1e-5 # iron (should really be 2.5e-1)
else:
values[0] = -6.4e-6 # copper (yes, it's negative!)
mu = Permeability(mesh, degree=1)
# Define variational problem
A_z = TrialFunction(V)
v = TestFunction(V)
a = (1 / mu)*dot(grad(A_z), grad(v))*dx
L_N = sum(J_N*v*dx(i) for i in range(2, 2 + n))
L_S = sum(J_S*v*dx(i) for i in range(2 + n, 2 + 2*n))
L = L_N + L_S
# Solve variational problem
A_z = Function(V)
solve(a == L, A_z, bc)
# Compute magnetic field (B = curl A)
W = VectorFunctionSpace(mesh, 'P', 1)
B = project(as_vector((A_z.dx(1), -A_z.dx(0))), W)
# Plot solution
plot(A_z)
plot(B)
# Save solution to file
vtkfile_A_z = File('magnetostatics/potential.pvd')
vtkfile_B = File('magnetostatics/field.pvd')
vtkfile_A_z << A_z
vtkfile_B << B
# Hold plot
interactive()