Modeling Voltage-driven Coils

NGSolve
1

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 × [b]a[/b] , curl × [b]a[/b]′)_Ω+ [i]I[/i] ( [b]t[/b], [b]a[/b]′)_Ω=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 [b]a[/b] \cdot [b]t[/b] + R [i]I[/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.

2

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”)