Dear community,
I’m trying to compute the eddy currents inside a conductive plate, with an imposed current inside a coil. After reading quite a few references on this topic, I ended up modeling the time-harmonic Maxwell equations with the magnetic vector potential A and the electric potential \varphi defined by \nabla \times A = B and E=j\omega \sigma (A + \nabla \varphi). To illustrate my problem, I’ve reduced the geometry to a unit cube, where both the coil and the plates are boxes with boundaries conforming to the structured mesh, and set all physical parameters to 1, except for the conductivity that is non-null on the plate and 0 elsewhere.
The variational forms for the eddy-currents problem is then the following :
The last two terms in a represent the Coulomb gauge \nabla \cdot A = 0 written as a saddle point problem, required to have a well-posed problem (see here). I discretize this problem with 1st kind Nedelec elements for A and first order Lagrange elements for \varphi, and use homogeneous boundary conditions n \times A = 0 and \varphi = 0 on \partial \Omega. As an example. the source currents is simply a constant unit vector along the x axis. I am then computing the eddy currents from the solution as J_\text{eddy} = j\sigma(A+\nabla \varphi) inside the plate and plot it to see the solution.
Here’s my code for this simplified problem:
import numpy as np
from mpi4py import MPI
import dolfinx, dolfinx.io, dolfinx.fem as fem, dolfinx.mesh
from petsc4py import PETSc
import ufl
cell_type = dolfinx.mesh.CellType.hexahedron
mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, 20, 20, 20, cell_type=cell_type)
mesh.topology.create_connectivity(2, 3)
# discretization space
if cell_type == dolfinx.mesh.CellType.hexahedron:
nedelec_string = "NCE"
if cell_type == dolfinx.mesh.CellType.tetrahedron:
nedelec_string = "N1curl"
Hcurl = ufl.FiniteElement(nedelec_string, mesh.ufl_cell(), degree=1)
H1 = ufl.FiniteElement("CG", mesh.ufl_cell(), degree=1)
W = fem.FunctionSpace(mesh, ufl.MixedElement([Hcurl, H1]))
N1curl = W.sub(0).collapse()[0]
CG1 = W.sub(1).collapse()[0]
A, phi = ufl.TrialFunctions(W)
N, psi = ufl.TestFunctions(W)
# define subdomains and meshteags for coil and plate
def coil_subdomain(x):
return ((0.45 <= x[1]) & (x[1] <= 0.56)) & ((0.35 <= x[2]) & (x[2] <= 0.41)) & ((0.3 <= x[0]) & (x[0] <= 0.7))
def plate_subdomain(x):
return ((0.1 <= x[1]) & (x[1] <= 0.91)) & ((0.2 <= x[2]) & (x[2] <= 0.301)) & ((0.1 <= x[0]) & (x[0] <= 0.91))
coil_entities = dolfinx.mesh.locate_entities(mesh, mesh.topology.dim, coil_subdomain)
plate_entities = dolfinx.mesh.locate_entities(mesh, mesh.topology.dim, plate_subdomain)
mt = dolfinx.mesh.meshtags(
mesh,
mesh.topology.dim,
np.concatenate([coil_entities, plate_entities]),
np.concatenate([np.full_like(coil_entities, 1), np.full_like(plate_entities, 2)]),
)
# export to plot problem geometry
with dolfinx.io.XDMFFile(mesh.comm, "mt.xdmf", "w") as file:
file.write_mesh(mesh)
file.write_meshtags(mt)
J = fem.Constant(mesh, np.array([1., 0, 0], dtype=PETSc.ScalarType))
# define conductivity in the plate only
DG0_scalar = fem.FunctionSpace(mesh, ("DG", 0))
sigma = fem.Function(DG0_scalar)
sigma.vector.set(0.)
sigma.x.array[plate_entities] = np.full(len(plate_entities), 1e5)
# variationnal problem
dx = ufl.Measure("dx", subdomain_data=mt)
a = ufl.inner(ufl.curl(A), ufl.curl(N)) * ufl.dx \
+ ufl.inner(A, ufl.grad(psi)) * ufl.dx \
+ ufl.inner(ufl.grad(phi), N) * ufl.dx \
+ 1j * sigma * ufl.inner(A + ufl.grad(phi), N + ufl.grad(psi)) * dx(2)
L = ufl.inner(J, N) * dx(1)
# dirichlet BCs
w_D = fem.Function(W)
def A_D_expr(x):
val = np.zeros_like(x, dtype=PETSc.ScalarType)
val[:, :] = 0.0
return val
def phi_D(x):
return np.full_like(x[0], 0.0)
w_D.sub(0).interpolate(A_D_expr)
w_D.sub(1).interpolate(phi_D)
boundary_facets = np.flatnonzero(dolfinx.mesh.compute_boundary_facets(mesh.topology))
boundary_dofs = fem.locate_dofs_topological(W.sub(1), mesh.topology.dim - 1, boundary_facets)
bc1 = fem.dirichletbc(w_D.sub(1), boundary_dofs)
boundary_dofs = fem.locate_dofs_topological(W.sub(0), mesh.topology.dim - 1, boundary_facets)
bc2 = fem.dirichletbc(w_D.sub(0), boundary_dofs)
# solve problem using a direct MUMPS solver
opts = {"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps" }
problem = fem.petsc.LinearProblem(a, L, bcs=[bc1, bc2], petsc_options=opts)
w_sol = problem.solve()
print(f"w_sol norm = {w_sol.vector.norm()}")
print(f"convergence reason {problem.solver.getConvergedReason()}")
# compute eddy currents and export solution
A_sol, phi_sol = w_sol.split()
DG1_vec = fem.VectorFunctionSpace(mesh, ("DG", 1), dim=3)
A_dg1 = fem.Function(DG1_vec, name="A")
A_dg1.interpolate(A_sol)
im = fem.Constant(mesh, 1j)
J_dg1 = fem.Function(DG1_vec, name="J")
J_dg1.interpolate(fem.Expression(im * sigma*(A_sol + ufl.grad(phi_sol)), DG1_vec.element.interpolation_points))
with dolfinx.io.VTXWriter(MPI.COMM_WORLD, "solution.bp", [A_dg1, J_dg1]) as f:
f.write(0.0)
The resolution of the linear system with the MUMPS direct solver works as expected, but when I plot the solution I observe a weird behaviour :
On the first plot, where I’m using a hexahedron based unit-cube, the computed eddy currents seem correct, however, on the second plot, where I’m using a tetrahedron based mesh, the solution seems wrong : on the surface of the plate the currents have non-null normal component which is non-physical (currents can’t go out of the plate). What’s weird is that I observe different solutions with different cells type for the same problem.
My guess would be that I am missing some interface condition that I need to enforce, but I can’t manage to figure out how to properly do it. Still, I find weird that I observe such differences when changing the mesh, as my problem is supposed to be well posed.
If anyone can provide me any hint that could point me towards what I’m doing wrong, it would be very helpful.
Thanks
