1. Problem 7 - Asymmetrical Conductor with a Hole

https://www.compumag.org/wp/wp-content/uploads/2018/06/problem7.pdf

Solve Eddy-current problem

\[\begin{align*} \nabla\times(\nu\nabla\times\mathbf{A}) + j\omega\sigma \mathbf{A} = \mathbf{J} \end{align*}\]

with

\[\begin{align*} \mathbf{A}\times\mathbf{n} = \mathbf{0} \end{align*}\]

on the far boundary. The magnetic permability \(\mu_0 = 4\pi \cdot 10^{-7}\) H/m of vaccuum, and the electric conductivity of the aluminium plate is \(\sigma_{Al} = 35.26\cdot10^6\) S/m.

from netgen.occ import *
from ngsolve import *
from ngsolve.webgui import Draw
import matplotlib.pyplot as plt
%matplotlib widget
import numpy as np

1.1. Generate Geometry

1.1.1. The aluminum plate:

wpplate = WorkPlane(Axes((0,0,0),Z,X))
f = wpplate.Rectangle(0.294,0.294).MoveTo(0.018,0.018).Rectangle(0.108,0.108).Reverse().Face()
plate = f.Extrude(0.019)

plate.faces.maxh=0.04
plate.solids.name="plate"
plate.faces.col=(0.3,0.3,0.3,1)
Draw (plate);

1.1.2. The coil:

wp = WorkPlane(Axes(p=(0.294,0.000,0.049), n=Z, h=Y)) # origin in 

coil2di = wp.MoveTo(0.100,0.100).RectangleC(0.100,0.100).Offset(0.025).Face()
coil2di.edges.name="coili"
coil2do = wp.MoveTo(0.100,0.100).RectangleC(0.100,0.100).Offset(0.050).Face()
coil2d = coil2do-coil2di

cutout = wp.MoveTo(0.100,0.100).RectangleC(0.300,0.100).Face() + \
         wp.MoveTo(0.100,0.100).RectangleC(0.100,0.300).Face()

coil2d = Glue( [coil2d*cutout, coil2d-cutout])

# Draw (coil2d)
coil = coil2d.Extrude(0.100)
coil.solids.name="coil"
Draw (coil);

airbox = Box(Pnt(-.2, -.2, -.2), Pnt(0.5, 0.5, 0.5))
airbox.faces.name="outer"
airbox.solids.name="air"

airbox.faces.col = (0,0,1,0.3)

airbox = airbox-coil-plate
shape = Glue([ coil, plate, airbox])
Draw (shape);

geo = OCCGeometry(shape)

xi_msm = np.array([0, 18, 36, 54, 72, 90, 108, 126, 144, 162, 180, 198, 216, 234, 252, 270, 288]) * 1e-3
xi_sim = np.linspace(0,0.288,500)

mesh = Mesh(geo.GenerateMesh(maxh=0.2))
mesh.Curve(5)
print("mats", set(mesh.GetMaterials()))
print("number of elements:", mesh.ne)

clipping_settings={"Clipping":{"enable":True, "z":-1, "dist":0}}
Draw(mesh, settings = clipping_settings);

mats {'air', 'plate', 'coil'}
number of elements: 14678

1.2. Coil currents

coilrect_xmin, coilrect_xmax = 0.294 - 0.150, 0.294-0.050
coilrect_ymin, coilrect_ymax = 0.050, 0.150

def Project(val, minval, maxval):
    return IfPos(val-minval, IfPos(val-maxval, maxval, val), minval)

projx = Project(x, coilrect_xmin, coilrect_xmax)
projy = Project(y, coilrect_ymin, coilrect_ymax)
tau_coil = CF( (projy-y, x-projx, 0) )
tau_coil /= Norm(tau_coil)
pot_coil = CF( (0, 0, sqrt((projy-y)**2+(projx-x)**2)-0.05) )

tau_coil_only = mesh.MaterialCF({"coil.*":tau_coil}, default=CF((0, 0, 0)))

clipping_settings={"Clipping":{"enable":True, "z":-1, "dist":-0.05}, "Objects":{"Clipping Plane":False, "Vectors":True}, "Vectors":{"grid_size":100}}
clipping_settings.update({"camera":{"transformations":[{"type":"move", "dir":(0,0,1), "dist":1.9}]}})
Draw(tau_coil_only, mesh,settings=clipping_settings, draw_surf=False); # top right, mesh, settings=clipping_settings)

1.3. FEM

mu = 4*pi*1e-7
sigma = mesh.MaterialCF({"plate":35.26e6}, default=1)

turns = 2742

def Solve(f):
    omega = 2*pi*f
    fes = HCurl(mesh, order=5, complex=True, dirichlet="outer", gradientdomains="plate")
    print("free dofs", sum(fes.FreeDofs()))
    u, v = fes.TnT()
    
    a = BilinearForm(fes, symmetric=True, condense=True)
    a += 1/mu*curl(u)*curl(v)*dx
    a += 1j*omega*sigma * u*v*dx

    # volume current density
    # f = LinearForm(-turns/(0.025 * 0.100)*tau_coil*v * dx("coil.*", bonus_intorder=4))

    # boundary current density + div-free correction
    f = LinearForm(
            -turns/(0.100)*tau_coil*v.Trace() * ds("coili.*", bonus_intorder=4) \
            +turns/(0.025 * 0.100)*pot_coil*curl(v) * dx("coil.*", bonus_intorder=4))
    
    A = GridFunction(fes)
    pre = Preconditioner(a, type="bddc", inverse="sparsecholesky")
    # a.Assemble()
    # f.Assemble()
    # inv = solvers.CGSolver(mat=a.mat, pre=pre, printrates='\r', maxiter=200)
    # A.vec[:] = inv*f.vec
    solvers.BVP(bf=a, lf=f, gf=A, pre=pre, \
                solver=solvers.CGSolver, solver_flags={"plotrates": True, "tol" : 1e-12})
    
    B = curl(A)
    J = -1j*omega*sigma * A
    return {"A":A, "B":B, "J":J}

with TaskManager():
    ret = Solve(f=50)
B, J = ret["B"], ret["J"]

1.4. Draw the B and J field

# clipping_settings={"Clipping":{"enable":True, "y":1, "z":0, "dist":0.1}, "Objects":{"Clipping Plane":False, "Vectors":True}, "Vectors":{"grid_size":100}}
# clipping_settings.update({"camera":{"transformations":[{"type":"rotateX", "angle":-90}, {"type":"move", "dir":(0,0,1.5), "dist":1}]}})
# Draw(B.real, mesh, settings=clipping_settings, max = 10e-3, min = 0, draw_surf=False);

clipping = { "function" : False,  "pnt" : (0,0.1,0), "vec" : (0,1,0) }
vectors = {"grid_size" : 50, "offset" : 0.5 }
Draw(B.real, mesh, clipping=clipping, vectors=vectors, max = 10e-3, min = 0, draw_surf=False);

from ngsolve.webgui import FieldLines
fieldlines = FieldLines (B.real, mesh.Materials(".*"), length=0.6, num_lines=100)
Draw (B.real, mesh, objects=[fieldlines], settings={"Objects":{"Surface":False, "Wireframe":False}}, \
     euler_angles=[-60,4,38]);

# clipping_settings={"Clipping":{"enable":True, "y":0, "z":-1, "dist":-0.132}, "Objects":{"Clipping Plane":False, "Vectors":True}, "Vectors":{"grid_size":150}}
# clipping_settings.update({"camera":{"transformations":[{"type":"move", "dir":(0,0,1), "dist":2.2}]}})
# Draw(mesh.MaterialCF({"plate":J.imag}, default = CF((0, 0, 0))), mesh, settings = clipping_settings, min = 0, max = 8e5);

clipping = { "function" : False,  "pnt" : (0,0,0.0189), "vec" : (0,0,-1) }
vectors = {"grid_size" : 100 }
Draw(mesh.MaterialCF({"plate":J.imag}, default = CF((0, 0, 0))), 
     mesh, clipping=clipping, vectors=vectors, min = 0, max = 8e5);

1.5. Reference values of B and J from the bench-mark

Bz_A1B1_50Hz_ref = np.array([ -4.9  -1.16j, -17.88 +2.48j, -22.13 +4.15j, -20.19 +4.j,   -15.67 +3.07j,
   0.36 +2.31j,  43.64 +1.89j,  78.11 +4.97j,  71.55+12.61j , 60.44+14.15j,
  53.91+13.04j,  52.62+12.4j,  53.81+12.05j, 56.91+12.27j,  59.24+12.66j,
  52.78 +9.96j,  27.61 +2.26j])
               
Bz_A2B2_50Hz_ref = np.array([-1.83-1.63j, -8.5-0.6j, -13.6-0.43j, -15.21+0.11j, -14.48+1.26j, -5.62+3.4j,
 28.77+6.53j, 60.34+10.25j, 61.84+11.83j, 56.64+11.83j, 53.4+11.01j, 52.36+10.58j, 53.93+10.8j, 56.82+10.54j, 
 59.48+10.62j, 52.08+9.03j, 26.56+1.79j])

Jy_A3B3_50Hz_ref = np.array([0.249-0.629j,  0.685-0.873j,  0+0j,  0+0j,  0+0j,  0+0j,  0+0j,  -0.015-0.593j,
             -0.103-0.249j,  -0.061-0.101j,  -0.004-0.001j,  0.051+0.087j,  0.095+0.182j,  0.135+0.322j,  
             0.104+0.555j,  -0.321+0.822j,  -0.687+0.855j,])


Jy_A4B4_50Hz_ref = np.array([0.461-0.662j,  0.621-0.664j,  0+0j,  0+0j,  0+0j,  0+0j,  0+0j,  1.573-1.027j,
             0.556-0.757j,  0.237-0.364j,  0.097-0.149j,  -0.034+0.015j,  -0.157+0.154j,  -0.305+0.311j,
             -0.478+0.508j,  -0.66+0.747j,  -1.217+1.034j])

1.5.1. Bz on A1 to B1

Bz_A1B1_sim = np.array([B[2](mesh(x, 0.072, 0.034)) for x in xi_sim])

plt.figure(1)
plt.clf()
plt.title("Evaluate B.z on Line A1 - B1")
plt.plot(xi_sim, -1e4*Bz_A1B1_sim.real, "b", label="sim Bz.real")
plt.plot(xi_sim, 1e4*Bz_A1B1_sim.imag, "r", label="sim Bz.imag")
plt.plot(xi_msm, Bz_A1B1_50Hz_ref.real, "x", label="ref Bz.real")
plt.plot(xi_msm, Bz_A1B1_50Hz_ref.imag, "x", label="ref Bz.imag")
plt.ylabel("B in Gauss = $10^{-4}$ T")
plt.legend()
plt.show()

1.5.2. Bz on A2 to B2

Bz_A2B2_sim = np.array([B[2](mesh(x, 0.144, 0.034)) for x in xi_sim])

plt.figure(2)
plt.clf()
plt.title("Evaluate B.z on Line A2 - B2")
plt.plot(xi_sim, -1e4 * Bz_A2B2_sim.real, "b", label="sim Bz.real")
plt.plot(xi_sim, 1e4 * Bz_A2B2_sim.imag, "r", label="sim Bz.imag")
plt.plot(xi_msm, Bz_A2B2_50Hz_ref.real, "b--x", label="ref Bz.real")
plt.plot(xi_msm, Bz_A2B2_50Hz_ref.imag, "r--x", label="ref Bz.imag")
plt.ylabel("B in Gauss = $10^{-4}$ T")
plt.legend()
plt.show()

clipping = { "function" : True,  "pnt" : (0,0.072,0), "vec" : (0,1,0) }

Draw (J[1].real, mesh, clipping = clipping, min=-0.5e6, max=0.5e6 )
Draw (J[1].imag, mesh, clipping = clipping, min=-0.5e6, max=0.5e6 );

Jy_A3B3_sim = np.array([1e-6*1j*J[1](mesh(x, 0.072, 0.019-1e-5)) for x in xi_sim])
Jy_A4B4_sim = np.array([1e-6*1j*J[1](mesh(x, 0.072, 0.000+1e-5)) for x in xi_sim])

1.6. Jy on A3 to B3 - plate top:

plt.figure(3)
plt.clf()
plt.title("Evaluate J.y on Line A3 - B3")
plt.plot(xi_sim, Jy_A3B3_sim.real, "b", label="sim Jy.real")
plt.plot(xi_sim, Jy_A3B3_sim.imag, "r", label="sim Jy.imag")
plt.plot(xi_msm, Jy_A4B4_50Hz_ref.real, "bx", label="ref Jy.real")
plt.plot(xi_msm, Jy_A4B4_50Hz_ref.imag, "rx", label="ref Jy.imag")
plt.ylabel("J in $10^{6}$ A/mm$^2$")
plt.ylim([-1.5,1.5])
plt.legend()
plt.show()

1.7. Jy on A4 to B4 - plate bottom

plt.figure(4)
plt.clf()
plt.title("Evaluate J.y on Line A4 - B4")
plt.plot(xi_sim, Jy_A4B4_sim.real, "b", label="sim Jy.real")
plt.plot(xi_sim, Jy_A4B4_sim.imag, "r", label="sim Jy.imag")
plt.plot(xi_msm, Jy_A3B3_50Hz_ref.real, "bx", label="ref Jy.real")
plt.plot(xi_msm, Jy_A3B3_50Hz_ref.imag, "rx", label="ref Jy.imag")
plt.ylabel("J in $10^{6}$ A/mm$^2$")
plt.ylim([-1.2,1.2])
plt.legend()
plt.show()