Hi guys!
I’m trying to solve a very simple PDE with anisotropic material property: The PDE is fairly normal magnetostatic equation that solves for magnetic vector potential
\nabla \times \upsilon \nabla \times \bf{A} = \bf{J}
where the reluctivity matrix can be defined as inverse of permeability matrix:
\upsilon = \mu ^{-1} =
\begin{bmatrix}
\upsilon_{x} & 0 & 0 \\
0 & \upsilon_{y} & 0 \\
0 & 0 & \upsilon_{z}
\end{bmatrix}
=
\begin{bmatrix}
{1}/{\mu_{x}} & 0 & 0 \\
0 & {1}/{\mu_{y}} & 0 \\
0 & 0 & {1}/{\mu_{z}}
\end{bmatrix}
As for electromagnetism, the higher the magnetic permeability (lower reluctivity) the higher the magnetic should be. In this case, where there is such high anisotropy, the flux should align itself along X axis.
But the thing is, I’m getting totally opposite effect where flux aligns itself along Y axis!
I’m not sure where does this effect come from… I’ve simplified example to 2D problem.
Am I doing something wrong here or I do not understand how FEniCS operates on gradients, tensors, coordinates?
from dolfin import *
from mshr import *
from scipy import constants
import matplotlib.pyplot as plt
domain = Rectangle(Point(-10, -10), Point(10, 10))
mesh = generate_mesh(domain, 64)
plot(mesh)
plt.show()
# function space
V = FunctionSpace(mesh, 'P', 1)
# boundary conditions
walls = 'on_boundary && (near(abs(x[0]), 10.0) || near(abs(x[1]), 10.0))'
bc = DirichletBC(V, Constant(0.0), walls)
tol = 1e-14
# Wire
class Omega_0(SubDomain):
def inside(self, x, on_boundary):
return x[0]**2 + x[1]**2 <= 4 + tol
# Space
class Omega_1(SubDomain):
def inside(self, x, on_boundary):
return x[0]**2 + x[1]**2 > 4 + tol
materials = MeshFunction('size_t', mesh, mesh.topology().dim())
subdomain_0 = Omega_0()
subdomain_1 = Omega_1()
subdomain_0.mark(materials, 0)
subdomain_1.mark(materials, 1)
plot(materials)
plt.show()
dx = Measure('dx', domain = mesh, subdomain_data = materials)
A_z = Function(V) # magnetic vector potential
v = TestFunction(V)
J = 5e6; # current source
# anisotropic material parameters
reluctivity = as_matrix( ( \
( 1/(constants.mu_0*10000), 0 ), \
( 0, 1/(constants.mu_0*1) ) \
) )
"""
reluctivity = 1/constants.mu_0
"""
F = inner(reluctivity*grad(A_z), grad(v))*dx - J*v*dx(0)
solve(F == 0, A_z, bc)
plot(A_z)
plt.xlabel('x')
plt.ylabel('y')
plt.show()
W = VectorFunctionSpace(mesh, 'P', 1)
Bx = A_z.dx(1)
By = -A_z.dx(0)
B = project( as_vector(( Bx, By )), W )
plot(B)
plt.xlabel('x')
plt.ylabel('y')
plt.show()