This page was generated from wta/coil.ipynb.
Magnetostatics¶
[1]:
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:¶
[2]:
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
[3]:
DrawGeo (coil);
[4]:
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:¶
[5]:
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:
[6]:
mesh.ne, mesh.nv, mesh.GetMaterials(), mesh.GetBoundaries()
[6]:
(102257,
17615,
('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
[7]:
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
[8]:
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\]
[9]:
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: 530927
[10]:
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.01096974634896
CG iteration 2, residual = 0.07848014588046139
CG iteration 3, residual = 0.014970621562618656
CG iteration 4, residual = 0.007249295251503955
CG iteration 5, residual = 0.0037126847844791686
CG iteration 6, residual = 0.002138972780124001
CG iteration 7, residual = 0.0012849799423056958
CG iteration 8, residual = 0.0008978901627578391
CG iteration 9, residual = 0.0006987952873078662
CG iteration 10, residual = 0.000490804689818925
CG iteration 11, residual = 0.00035040084437774205
CG iteration 12, residual = 0.0002809206411471045
CG iteration 13, residual = 0.00022174582794228542
CG iteration 14, residual = 0.00016575592975271304
CG iteration 15, residual = 0.00012076112480096311
CG iteration 16, residual = 8.602456203835686e-05
CG iteration 17, residual = 5.77383319791224e-05
CG iteration 18, residual = 3.641976854236137e-05
CG iteration 19, residual = 2.256718543934991e-05
CG iteration 20, residual = 1.381955166754894e-05
CG iteration 21, residual = 9.419863586296443e-06
CG iteration 22, residual = 6.349325209121514e-06
CG iteration 23, residual = 5.062514923721583e-06
CG iteration 24, residual = 2.8663882557314465e-06
CG iteration 25, residual = 1.7565387857869199e-06
CG iteration 26, residual = 1.0930621419706796e-06
CG iteration 27, residual = 6.781840691442851e-07
CG iteration 28, residual = 4.504279787363252e-07
CG iteration 29, residual = 2.9311353143049714e-07
CG iteration 30, residual = 1.946186147057951e-07
CG iteration 31, residual = 1.3460684082528044e-07
CG iteration 32, residual = 1.0170724331880104e-07
CG iteration 33, residual = 7.996339634571439e-08
CG iteration 34, residual = 7.390847592792018e-08
CG iteration 35, residual = 4.700535231975139e-08
CG iteration 36, residual = 3.3959799538322095e-08
CG iteration 37, residual = 2.443058182611789e-08
CG iteration 38, residual = 1.7477871591928628e-08
CG iteration 39, residual = 1.263929866463822e-08
CG iteration 40, residual = 8.418998837120271e-09
CG iteration 41, residual = 5.493013325268904e-09
CG iteration 42, residual = 3.5183817755220558e-09
CG iteration 43, residual = 2.348712558171994e-09
CG iteration 44, residual = 1.9249775711465884e-09
CG iteration 45, residual = 1.36960225023062e-09
CG iteration 46, residual = 9.056441135977306e-10
CG iteration 47, residual = 6.300064706306453e-10
CG iteration 48, residual = 4.097890091014912e-10
CG iteration 49, residual = 2.7184061051551916e-10
CG iteration 50, residual = 1.7825961644221675e-10
CG iteration 51, residual = 1.2239946630406048e-10
CG iteration 52, residual = 8.519824152999747e-11
CG iteration 53, residual = 6.238229136141398e-11
CG iteration 54, residual = 4.867286960008958e-11
CG iteration 55, residual = 4.3259563963704654e-11
CG iteration 56, residual = 2.82494903170851e-11
CG iteration 57, residual = 2.10638326397284e-11
[11]:
Draw (curl(gfu), mesh, draw_surf=False, \
min=0, max=3e-4, clipping = { "y":1, "z" : 0, "function":False}, vectors = { "grid_size":100});
[ ]: