Transient electromagnetic PDE problem to calculate eddy currents

variational formulation
1

I am trying to create a script in FEniCS that will calculate electromagnetic fields with eddy currents taken into consideration. The domain was divided into two subdomains: copper wire (yellow), and air (purple). For simplicity both have the same permeability, but only copper has conductivity of \sigma = 58~\textrm{MS/m}. The conductivity of air is not taken into account, but i know that for some cases it should be taken into consideration.
The boundary condition is \textbf{A} = 0.

Domain

The mesh for this problem is shown below:

Mesh

The governing PDE for transient electromagnetic problem with eddy currents I have decided to use is based on magnetic vector potential \bf{A}

\nabla \times \frac{1}{\mu_0} \nabla \times \textbf{A} = -\sigma \frac{\partial \textbf{A}}{\partial t} + \textbf{J}_s(t)

where: \textbf{J}_s is the source current and term \sigma \frac{\partial \textbf{A}}{\partial t} is the eddy current \textbf{J}_e.

The partial derivative over time is represented by finite difference: \frac{\textbf{A}_n - \textbf{A}_{n-1}}{\Delta t}, where \textbf{A}_n is distribution of magnetic vector potential at current step and \textbf{A}_{n-1} at previous step. The \Delta t is the time step.

This gives the PDE as:
\nabla \times \frac{1}{\mu_0} \nabla \times \textbf{A}_n + \sigma \frac{\textbf{A}_n}{\Delta t} = \sigma \frac{\textbf{A}_{n-1}}{\Delta t} + \textbf{J}_s(t)

and weak formulation with test function \bf{v}:
\int_{\Omega} (\frac{1}{\mu_0} \nabla \times \textbf{A}_n) \cdot (\nabla \times \textbf{v}) d\Omega + \int_{\Omega_{Cu}} \sigma \frac{\textbf{A}_n}{\Delta t} \textbf{v} d\Omega_{Cu} = \int_{\Omega_{Cu}} \sigma \frac{\textbf{A}_{n-1}}{\Delta t} \textbf{v} d\Omega_{Cu} + \int_{\Omega_{Cu}} \textbf{J}_s(t) \textbf{v} d\Omega_{Cu}

What I am expecting to get is some sort of current distribution that relates to skin depth of copper at frequency of 50 kHz, but this doesn’t seem to happen and I don’t know why.

In loop I’m updating the current density, which follows a sinusoid, and assigning field distribution to be used as the one calculated in previous.
The time derivative of \bf{A} is calculated only in region of copper (yellow). The eddy currents \bf{J}_e are calculated separately from source currents \bf{J}_s, and combined into total current density \bf{J}_{tot} to be saved in .pvd file.

Could someone please point out how to fix this issue, assemble equations or find mistakes? At frequency of 50 kHz, the skin effect should be quite visible at given diameter of circle, yet it isn’t.

from dolfin import *
import mshr
import matplotlib.pyplot as plt
import numpy as np

pvd_A_z = File('A_z.pvd')
pvd_B_vec = File('B_vec.pvd')
pvd_B_mag = File('B_mag.pvd')
pvd_J_e = File('J_e.pvd')
pvd_J_tot = File('J_tot.pvd')

def curl2D(u):
	return as_vector( (u.dx(1), -u.dx(0)) )
	
# define constants
mm = 1e-3
mu = 4*np.pi*1e-7 # H/m
sigma_Cu = 58e6 # MS/m

f = 50e3 # Hz
T = 1/f # s
nSteps = 20
dt = T/nSteps

# create geometry
region = mshr.Rectangle(Point(-10*mm,-10*mm), Point(10*mm,10*mm))
wire = mshr.Circle(Point(0*mm,0*mm), 3*mm)

region.set_subdomain(1, wire)

# mesh and subdomains
mesh = mshr.generate_mesh(region, 64)
markers = MeshFunction('size_t', mesh, 2, mesh.domains())

dx = Measure('dx', domain=mesh, subdomain_data=markers)

# function space
V = FunctionSpace(mesh, 'P', 1)
W = VectorFunctionSpace(mesh, 'P', 1)

# boundary conditions
bc = DirichletBC(V, Constant(0), 'on_boundary')

# Define problem
A_n = TrialFunction(V) # current step
A_n_1 = Function(V) # previous step
v = TestFunction(V)
J = Expression('Jmax', Jmax=0.0, degree=0) # source

a = inner((1/mu)*curl2D(A_n), curl2D(v))*dx \
	+ (sigma_Cu/dt)*A_n*v*dx(1)
L = J*v*dx(1) + (sigma_Cu/dt)*A_n_1*v*dx(1)

A_n = Function(V) # current step

delA = Function(V) # Difference between A_n and A_n_1
J_tot = Function(V) # total current density

# Main loop
t = 0.0
for i in range(nSteps):
	solve(a == L, A_n, bc)
	
	A_n.rename("A_z", "1")
	pvd_A_z << (A_n, t)
	
	# Project source currents
	J_s = TrialFunction(V)
	w = TestFunction(V)
	a_h = J_s*w*dx
	L_h = J*w*dx(1)
	J_s = Function(V)
	solve(a_h == L_h, J_s)
	
	# Calculate eddy currents
	J_e = TrialFunction(V)
	w = TestFunction(V)
	delA.vector()[:] = A_n.vector()[:] - A_n_1.vector()[:]
	a_h = J_e*w*dx
	L_h = (sigma_Cu/dt)*delA*w*dx(1)
	J_e = Function(V)
	solve(a_h == L_h, J_e)
	J_e.rename("J_e", "1")
	pvd_J_e << (J_e, t)
	
	# Combine currents
	J_tot.vector()[:] = J_e.vector()[:] + J_s.vector()[:]
	J_tot.rename("J_tot", "1")
	pvd_J_tot << (J_tot, t)
	
	# Get B field
	B = project(curl2D(A_n), W)
	B.rename("B_vec", "1")
	pvd_B_vec << (B, t)
	
	# Update
	t += dt
	J.Jmax = 2e6*sin(2*np.pi*f*t)
	A_n_1.assign(A_n)
2 Likes