Closed coil magnetostatic problem

NGSolve
1

Hi,

I am trying to follow the magnetostatic example in the I-tutorials but I want to close the coil in on itself and still solve for the current in the coil. (And the magnetic field around the coil)

Does anyone have an idea how I achieve this? In essence I have to define the dirichlet boundary condition and current source on either side of the same surface where the two ends meet.

The link to the I-tutorial for reference: Magnetostatics — NGS-Py 6.2.2604 documentation

Thanks in advance!

2

Hi,
a nice way is to use a mixed formulation, and search for the coil current in H(div). On the interface, you can put the integral \int U tau*n ds to the right hand side, where U is the applied voltage over the interface.
Joachim

3

Hi,

Thank you so much for the reply! I am a complete beginner when it comes to ngsolve but tried to do something like what you describe which I have appended below. I however get that the iterations just return nan (CG iteration 1, residual = nan, CG iteration 2, residual = nan, etc)

I am not sure exactly what it is that I am missing?

Thanks again,

Hampus

import ngsolve as ngs
import netgen.occ as occ

spine = occ.Wire(occ.Circle((0,0,0), (0,0,1), 1))
circle = occ.Face(occ.Wire([occ.Circle((1,0,0), (0,1,0), 0.1)]))
circle.faces.name = "cut"
coil = occ.Glue([occ.Pipe(spine, circle), circle])
coil.solids.name = "coil"
coil.faces.name = "outer"
geometry = occ.OCCGeometry(coil)

# Mesh
mesh = ngs.Mesh(geometry.GenerateMesh())

# FEM
U = 1       # Voltage
sigma = 10  # conductivity
divfes = ngs.HDiv(mesh, order=3, dirichlet="outer")
l2fes = ngs.L2(mesh, order=2)
fes = divfes*l2fes
(ud, ul2), (vd, vl2) = fes.TnT()
a = ngs.BilinearForm((1/sigma*ud*vd + ngs.div(ud)*vl2 + ngs.div(vd)*ul2)*ngs.dx)
c = ngs.Preconditioner(a, type="bddc")
n = ngs.specialcf.normal(mesh.dim)
f = ngs.LinearForm(U*(vd.Trace()*n*ngs.ds(f"cut")))

gfm = ngs.GridFunction(fes)
with ngs.TaskManager():
    a.Assemble()
    f.Assemble()
    ngs.solvers.CG(sol=gfm.vec, rhs=f.vec, mat=a.mat, pre=c.mat)
4

Here is a working code:

  • check your geometry, Draw it, look if boundary labels show up (sometimes tricky)
  • I used H(div) with penalty, no saddle point problem (this is next step, if you like)
  • start with direct solvers, preconditioners later
import ngsolve as ngs
import netgen.occ as occ
from ngsolve.webgui import Draw

spine = occ.Wire(occ.Circle((0,0,0), (0,0,1), 1))
circle = occ.Face(occ.Wire([occ.Circle((1,0,0), (0,1,0), 0.1)]))
circle.faces.name = "cut"
coil = occ.Glue([occ.Pipe(spine, circle), circle])
coil.solids.name = "coil"
for f in coil.faces:
    if f.name==None:
        f.name = "outer"
geometry = occ.OCCGeometry(coil)

# Mesh
mesh = ngs.Mesh(geometry.GenerateMesh()).Curve(3)
Draw (coil)
Draw (mesh);
print (mesh.GetBoundaries())



U = 1       # Voltage
sigma = 10  # conductivity
fes = ngs.HDiv(mesh, order=2, dirichlet="outer")
ud,vd = fes.TnT()
a = ngs.BilinearForm(1/sigma*ud*vd*ngs.dx + 1e6 * ngs.div(ud)*ngs.div(vd)*ngs.dx)
n = ngs.specialcf.normal(mesh.dim)
f = ngs.LinearForm(U*(vd.Trace()*n*ngs.ds("cut")))
gfm = ngs.GridFunction(fes)


with ngs.TaskManager():
    a.Assemble()
    f.Assemble()
    gfm.vec.data = a.mat.Inverse(freedofs=fes.FreeDofs())*f.vec

Draw (gfm, draw_surf=False, clipping = { "pnt" : (0,0,0), "vec" : (0,0,-1), "function" : True }, min=1, max=2);
Draw (gfm, draw_surf=False, clipping = { "pnt" : (0,0,0), "vec" : (0,0,-1) }, vectors = { "grid_size" : 100 }, min=1, max=2);
5

Thank you for the help! It works like a charm!

Is it also possible to somehow treat two sides of the same surface as different surfaces?

For instance if you wanted to solve the above problem in H1:

fes = ngs.H1(mesh, order=3, dirichlet="cut_right") # <-- one side of the cut
u, v = fes.TnT()
a = ngs.BilinearForm(sigma*grad(u)*grad(v)*ngs.dx)
f = ngs.LinearForm(J_in*v*ds("cut_left")) # <-- other side of the cut

Or if you were interested in stipulating different dirichlet conditions on either side?

Hampus

6

You can do a similar technique to here:

Best
Christopher