This page was generated from wta/coil.ipynb.

Magnetostatics

from netgen.occ import *
from ngsolve import *
from ngsolve.webgui import Draw
from netgen.webgui import Draw as DrawGeo
import math

model of the coil:

cyl = Cylinder((0,0,0), Z, r=0.01, h=0.03).faces[0]
heli = Edge(Segment((0,0), (12*math.pi, 0.03)), cyl)
ps = heli.start
vs = heli.start_tangent
pe = heli.end
ve = heli.end_tangent

e1 = Segment((0,0,-0.03), (0,0,-0.01))
c1 = BezierCurve( [(0,0,-0.01), (0,0,0), ps-vs, ps])
e2 = Segment((0,0,0.04), (0,0,0.06))
c2 = BezierCurve( [pe, pe+ve, (0,0,0.03), (0,0,0.04)])
spiral = Wire([e1, c1, heli, c2, e2])
circ = Face(Wire([Circle((0,0,-0.03), Z, 0.001)]))
coil = Pipe(spiral, circ)

coil.faces.maxh=0.2
coil.faces.name="coilbnd"
coil.faces.Max(Z).name="in"
coil.faces.Min(Z).name="out"
coil.mat("coil")
crosssection = coil.faces.Max(Z).mass

DrawGeo (coil);

box = Box((-0.04,-0.04,-0.03), (0.04,0.04,0.06))
box.faces.name = "outer"
air = box-coil
air.mat("air");

mesh-generation of coil and air-box:

geo = OCCGeometry(Glue([coil,air]))
with TaskManager():
    mesh = Mesh(geo.GenerateMesh(meshsize.coarse, maxh=0.01)).Curve(3)

Draw (mesh, clipping={"y":1, "z":0, "dist":0.012});

checking mesh data materials and boundaries:

mesh.ne, mesh.nv, mesh.GetMaterials(), mesh.GetBoundaries()

(154246,
 26567,
 ('coil', 'air'),
 ('out',
  'coilbnd',
  'coilbnd',
  'coilbnd',
  'coilbnd',
  'coilbnd',
  'in',
  'outer',
  'outer',
  'outer',
  'outer',
  'outer',
  'outer'))

Solve a potential problem to determine current density in wire:

on the domain \(\Omega_{\text{coil}}\): \begin{eqnarray*} j & = & \sigma \nabla \Phi \\ \operatorname{div} j & = & 0 \end{eqnarray*} port boundary conditions: \begin{eqnarray*} \Phi & = & 0 \qquad \qquad \text{on } \Gamma_{\text{out}}, \\ j_n & = & \frac{1}{|S|} \quad \qquad \text{on } \Gamma_{\text{in}}, \end{eqnarray*} and \(j_n=0\) else

fespot = H1(mesh, order=3, definedon="coil", dirichlet="out")
phi,psi = fespot.TnT()
sigma = 58.7e6
bfa = BilinearForm(sigma*grad(phi)*grad(psi)*dx).Assemble()
inv = bfa.mat.Inverse(freedofs=fespot.FreeDofs(), inverse="sparsecholesky")
lff = LinearForm(1/crosssection*psi*ds("in")).Assemble()
gfphi = GridFunction(fespot)
gfphi.vec.data = inv * lff.vec

Draw (gfphi, draw_vol=False, clipping={"y":1, "z":0, "dist":0.012});

Solve magnetostatic problem:

current source is current from potential equation:

\[\int \mu^{-1} \operatorname{curl} u \cdot \operatorname{curl} v \, dx = \int j \cdot v \, dx\]

fes = HCurl(mesh, order=2, nograds=True)
print ("HCurl dofs:", fes.ndof)
u,v = fes.TnT()
mu = 4*math.pi*1e-7
a = BilinearForm(1/mu*curl(u)*curl(v)*dx+1e-6/mu*u*v*dx)
pre = Preconditioner(a, "bddc")
f = LinearForm(sigma*grad(gfphi)*v*dx("coil"))
with TaskManager():
    a.Assemble()
    f.Assemble()

HCurl dofs: 799836

from ngsolve.krylovspace import CGSolver
inv = CGSolver(a.mat, pre, printrates=True)
gfu = GridFunction(fes)
with TaskManager():
    gfu.vec.data = inv * f.vec

CG iteration 1, residual = 23.416986574283737
CG iteration 2, residual = 0.06705058462253938
CG iteration 3, residual = 0.01734237793941035
CG iteration 4, residual = 0.011292648574682668
CG iteration 5, residual = 0.007228002610033505
CG iteration 6, residual = 0.004590707306041366
CG iteration 7, residual = 0.0030083060795949336
CG iteration 8, residual = 0.0021614820269268285
CG iteration 9, residual = 0.0013851987760443646
CG iteration 10, residual = 0.001017031052403044
CG iteration 11, residual = 0.0007711358816311045
CG iteration 12, residual = 0.0005753312522126517
CG iteration 13, residual = 0.00046120468765010065
CG iteration 14, residual = 0.00037992763898862385
CG iteration 15, residual = 0.00032197420958833697
CG iteration 16, residual = 0.0002727496356143563
CG iteration 17, residual = 0.00022256108308747258
CG iteration 18, residual = 0.00017269048957516996
CG iteration 19, residual = 0.0001305626294064451
CG iteration 20, residual = 0.00010072613836770141
CG iteration 21, residual = 7.28469164988659e-05
CG iteration 22, residual = 5.6186102669194734e-05
CG iteration 23, residual = 4.007243144790527e-05
CG iteration 24, residual = 2.914483874676209e-05
CG iteration 25, residual = 2.1388441353490006e-05
CG iteration 26, residual = 1.5207357282299595e-05
CG iteration 27, residual = 1.1164901036478488e-05
CG iteration 28, residual = 8.64375712310422e-06
CG iteration 29, residual = 6.782758325300258e-06
CG iteration 30, residual = 5.6964895398021706e-06
CG iteration 31, residual = 4.6404953318739245e-06
CG iteration 32, residual = 4.037425501745331e-06
CG iteration 33, residual = 3.548224384719152e-06
CG iteration 34, residual = 2.6939265273856113e-06
CG iteration 35, residual = 2.0778132505856147e-06
CG iteration 36, residual = 1.5896265408369998e-06
CG iteration 37, residual = 1.1986984393137937e-06
CG iteration 38, residual = 8.974825852510094e-07
CG iteration 39, residual = 6.570450504373847e-07
CG iteration 40, residual = 4.93251918653432e-07
CG iteration 41, residual = 3.7115560281778636e-07
CG iteration 42, residual = 3.2421641354038096e-07
CG iteration 43, residual = 2.441419052845636e-07
CG iteration 44, residual = 1.9403486252029256e-07
CG iteration 45, residual = 1.6576868634067552e-07
CG iteration 46, residual = 1.357905914960379e-07
CG iteration 47, residual = 1.0843013916602062e-07
CG iteration 48, residual = 8.307860190038368e-08
CG iteration 49, residual = 7.143203713933814e-08
CG iteration 50, residual = 5.5321032358664197e-08
CG iteration 51, residual = 4.1296563550275244e-08
CG iteration 52, residual = 3.065721991602448e-08
CG iteration 53, residual = 2.2885784297685955e-08
CG iteration 54, residual = 1.6850780371594167e-08
CG iteration 55, residual = 1.3679883407975996e-08
CG iteration 56, residual = 1.1039943134699962e-08
CG iteration 57, residual = 8.140490522387797e-09
CG iteration 58, residual = 6.00032663887363e-09
CG iteration 59, residual = 5.127430798212429e-09
CG iteration 60, residual = 4.319389077146737e-09
CG iteration 61, residual = 3.1745076721841292e-09
CG iteration 62, residual = 2.563182528250854e-09
CG iteration 63, residual = 2.1241876484500757e-09
CG iteration 64, residual = 1.6363791842597117e-09
CG iteration 65, residual = 1.2590517273802158e-09
CG iteration 66, residual = 9.504878471082639e-10
CG iteration 67, residual = 7.597383753594858e-10
CG iteration 68, residual = 5.970377567188569e-10
CG iteration 69, residual = 4.481707548413982e-10
CG iteration 70, residual = 3.45534506855563e-10
CG iteration 71, residual = 2.829156657820269e-10
CG iteration 72, residual = 2.181051946088125e-10
CG iteration 73, residual = 1.5982147570450668e-10
CG iteration 74, residual = 1.1959067138816952e-10
CG iteration 75, residual = 9.636718694975352e-11
CG iteration 76, residual = 7.340125327926563e-11
CG iteration 77, residual = 5.5080704340447666e-11
CG iteration 78, residual = 4.411192161647098e-11
CG iteration 79, residual = 3.510301620578237e-11
CG iteration 80, residual = 2.7577776929927507e-11
CG iteration 81, residual = 2.1845112147778327e-11

Draw (curl(gfu), mesh, draw_surf=False, \
      min=0, max=3e-4, clipping = { "y":1, "z" : 0, "function":False}, vectors = { "grid_size":100});