# Modeling Voltage-driven Coils

• philscher
• New Member
• 2 months 6 days ago - 2 months 6 days ago #4634
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.
Last edit: 2 months 6 days ago by philscher.

• philscher
• New Member
• 2 months 2 days ago - 2 months 2 days ago #4640
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

rect   = Rectangle(pmin=(-0.05,-0.05), pmax=(0.05,0.05), mat="Air", bc="Outer" )
air = rect - circle

geo = CSG2d()

from ngsolve import *

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 += -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 * I_ / CoilArea * dx("Coil")

S = GridFunction(X)

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.vec.data
n = n + 1
print(f"Time: {t} Current: {CoilCurrent}")
Current.extend(CoilCurrent)
time.extend([t])
Draw(S.components, 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")

Last edit: 2 months 2 days ago by philscher.