Forum Message

 

 

We have moved the forum to https://forum.ngsolve.org . This is an archived version of the topics until 05/05/23. All the topics were moved to the new forum and conversations can be continued there. This forum is just kept as legacy to not invalidate old links. If you want to continue a conversation just look for the topic in the new forum.

Notice

The forum is in read only mode.

Modeling Voltage-driven Coils

More
1 year 2 months ago - 1 year 2 months ago #4634 by philscher
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: 1 year 2 months ago by philscher.
More
1 year 2 months ago - 1 year 2 months ago #4640 by philscher
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")


 
Last edit: 1 year 2 months ago by philscher.
Time to create page: 0.133 seconds