Modeling Voltage-driven Coils

Dear NGSolve community,

 I wonder how to model a voltage driven (electromagnetic) stranded coil in ngsolve ?The weak-forms in the coil in the domain Ω are given by:

The field problem in the a-formulation

       ( nu curl × a , curl × a′)_Ω+ I ( t, a′)_Ω=0,

where a is the magnetic vector potential, I the (unknown) current flowing through each strand of the coil, and t is the electric current density in the coil when a current of 1A is applied.

The electric circuit problem (which introduces a single additional dof)

         ∫_Ω partial_t a \cdot t + R I = V,

where R is the specified resistance and V is the specified voltage.

Thanks a lot for any reply and many thanks for an amazing software.

I think we are able to bring in the additional dof through the NumberSpaces, e.g., 

V = H1(mesh, order = 1, dirichlet = "Outer")
N = NumberSpace(mesh)
X = FESpace([V,N])

I wonder however if there is a simpler way to setup the bilinear form for the circuit equation, which requires the division by the CoilArea:

a += tau/dt * A * I_ / CoilArea * dx("Coil") 

Any suggestions ? I have attached the file in case anybody is interested (or finds any improvements / bugs).

---

Adding file here as attaching is not working for my browser

from netgen.geom2d import CSG2d, Circle, Rectangle

CoilRadius = 0.01

circle = Circle( center=(0,0), radius=CoilRadius, mat="Coil", bc="Boundary" )
rect   = Rectangle(pmin=(-0.05,-0.05), pmax=(0.05,0.05), mat="Air", bc="Outer" )
air = rect - circle

geo = CSG2d()

geo.Add(circle)
geo.Add(air)

from ngsolve import *
mesh = Mesh(geo.GenerateMesh(maxh=2.e-2, quad_dominated=True))


nu = 1./1.257e-6

U = 5.
R = 2.
Turns = 20

dt = 5.e-2
t = 0
Tend = 2.

##################


V = H1(mesh, order = 1, dirichlet = "Outer")
N = NumberSpace(mesh)
X = FESpace([V,N])

(A,I), (A_,I_) = X.TnT()


CoilArea = Integrate(CoefficientFunction(1), mesh, definedon=mesh.Materials("Coil"))
tau = 1. / CoilArea * Turns

a = BilinearForm(X, symmetric=False)
a += nu*grad(A)*grad(A_) * dx 
a += -I * tau * A_ * dx("Coil")
a += tau/dt * A * I_ / CoilArea * dx("Coil") 
a += R * I * I_ / CoilArea * dx("Coil")

c = Preconditioner(a, "direct")

a.Assemble()

n = 0

Current = [ 0.]
time    = [ 0.]

Aold = GridFunction(X)

while t < Tend:

    f = LinearForm(X)
    f += U * I_ / CoilArea * dx("Coil")
    f += tau/dt * Aold.components[0] * I_ / CoilArea * dx("Coil") 

    S = GridFunction(X)

    with TaskManager():
            f.Assemble()
            bvp = BVP(bf=a, lf=f, gf=S, pre=c)
            bvp.Do()

    Aold.vec.data = S.vec.data
    t += dt

    CoilCurrent = S.components[1].vec.data
    n = n + 1
    print(f"Time: {t} Current: {CoilCurrent}")
    Current.extend(CoilCurrent)
    time.extend([t])
    Draw(S.components[0], mesh, "A") 


import pylab
import numpy

pylab.plot(time, Current)
pylab.plot(time, (U/R) * numpy.ones_like(time))

pylab.xlabel("Time")
pylab.ylabel("Current")
pylab.ylim((0, 3))

pylab.savefig("Current.png")