Force between magnet and coil (virtual work)

Domi
Topic Author
Offline
New Member

6 years 3 months ago - 6 years 3 months ago #1173 by Domi

Force between magnet and coil (virtual work) was created by Domi

Hello,
first, thank you for building this powerful software!
I am new to NGS and what I am trying to do is to compute the force between a permanent magnet and a rectangular coil. For that reason, I compute the energy in the system, displace the magnet and recompute the energy. The force should then result as F = (Energy1-Energy2)/displacement. But there is a problem with the calculation of the energy. The values are arbitrary huge.
Result: F = 1525527259.8982043, should be something like 5 N
I am grateful for any help!
Domi

Code:

from ngsolve import *
from netgen.csg import *
import netgen.gui

distance_magnet= 20.0/1000 #distance of the magnet for <20[mm] it is inside the coil
j=(1 * 700)/0.00005 #current density (I*N)/A [A/m^2]
step = 0.0001

w=[None, None]
box = OrthoBrick(Pnt(-0.090, -0.090, -0.090),Pnt(distance_magnet+0.105, 0.1246, 0.1084)).bc("outer")
coil1 = (OrthoBrick(Pnt(0, 0,0),Pnt(0.020, 0.0346, 0.0025)))
coil2 = (OrthoBrick(Pnt(0, 0, 0.0025),Pnt(0.020, 0.0025, 0.0159)))
coil3 = (OrthoBrick(Pnt(0, 0, 0.0159),Pnt(0.020, 0.0346, 0.0184)))
coil4 = (OrthoBrick(Pnt(0, 0.0321, 0.0025),Pnt(0.020, 0.0346, 0.0159)))

for y in range(0, 2):
    magnet = OrthoBrick(Pnt(distance_magnet, 0.0048, 0.0035),Pnt(distance_magnet + 0.015, 0.0298, 0.0145))
    air = box - magnet - coil1 - coil2 - coil3 - coil4
    geo = CSGeometry()
    geo.Add(air.mat("air"))
    geo.Add(coil1.mat("coil1").maxh(0.005))
    geo.Add(coil2.mat("coil2").maxh(0.005))
    geo.Add(coil3.mat("coil3").maxh(0.005))
    geo.Add(coil4.mat("coil4").maxh(0.005))
    geo.Add(magnet.mat("magnet").maxh(0.003))
    mesh = Mesh(geo.GenerateMesh(maxh=0.010))
    mesh.GetMaterials(), mesh.GetBoundaries()
    fes = HCurl(mesh, order=3, dirichlet="outer", flags = { "nograds" : True })
    print("ndof =", fes.ndof)
    u = fes.TrialFunction()
    v = fes.TestFunction()

    from math import pi
    mu0 = 4 * pi * 1e-7
    mur = CoefficientFunction([1.0445730167132 if mat == "magnet" else
                                1.0000000000000 if mat == "air" else 0.999991 #mu copper
                                for mat in mesh.GetMaterials()])

    a = BilinearForm(fes)
    a += SymbolicBFI(1 / (mu0 * mur) * curl(u) * curl(v) + 1e-6 / (mu0 * mur) * u * v)
    c = Preconditioner(a, "bddc")

    f = LinearForm(fes) 
    mag = CoefficientFunction((-837999.999999998,0,0)) #coercivity -> -837999.999999998 [A m^-1]
    f += SymbolicLFI(mag * curl(v), definedon=mesh.Materials("magnet"))
    f += SymbolicLFI(CoefficientFunction((0,j,0)) * v, definedon=mesh.Materials("coil1")) 
    f += SymbolicLFI(CoefficientFunction((0,0,-j)) * v, definedon=mesh.Materials("coil2"))
    f += SymbolicLFI(CoefficientFunction((0,-j,0)) * v, definedon=mesh.Materials("coil3"))
    f += SymbolicLFI(CoefficientFunction((0,0,j)) * v, definedon=mesh.Materials("coil4"))

    with TaskManager():
        a.Assemble()
        f.Assemble()

    gfu = GridFunction(fes)

    inv = CGSolver(a.mat, c.mat)
    gfu.vec.data = inv * f.vec
    Draw(curl(gfu), mesh, "B-field")

    w[y] = 0.5 * InnerProduct(gfu.vec, f.vec)
    print("W",y, "=" , w[y])
    distance_magnet += step

F = (w[0] - w[1])/step
print("F=", F)

Last edit: 6 years 3 months ago by Domi.

joachim
Offline
Administrator

6 years 3 months ago #1174 by joachim

Replied by joachim on topic Force between magnet and coil (virtual work)

Hi Domi,

just yesterday I stumbled over the same question:
The advanced method is not to modify the geometry and use finite differences, but using some calculus and compute the shape derivative.
It results in the so called Maxwell stress tensor, and its divergence ist the force density.
(One big problem with the FD is that two similar geometries may not result in two similar meshes.)

Some reference is here:
Franc?ois Henrotte, Kay Hameyer: Computation of electromagnetic force densities: Maxwell stress tensor vs. virtual work principle?
Journal of Computational and Applied Mathematics 168 (2004)

I have a DRAFT VERSION of the TEAM23 benchmark which shows how to compute the integral forces from the stress tensor. Not all terms are in yet and we don't get the reference values, but it shows how to do the force calculation in NGSolve with two different methods.

Best, Joachim

Attachments:

team23.py

The following user(s) said Thank You: creativeworker

Domi
Topic Author
Offline
New Member

6 years 2 months ago #1210 by Domi

Replied by Domi on topic Force between magnet and coil (virtual work)

Hey Joachim,
thank you for your fast and meaningful answer. Luckily my supervisor could help me understand it... Now it works great.
I compared the results to an ANSYS simulation and to several measurements. ANSYS delivers more or less exactly the same results. The forces measured are slightly lower but all in all also stunningly close to the calculated values.
The code is attached. There is a for-loop and an array "distances" to calculate the forces for different distances between coil and magnet.
Regards,
Domi

Attachments:

...Coil_magnet.py

Michi
New Member

2 years 2 months ago #4505 by Michi

Replied by Michi on topic Force between magnet and coil (virtual work)

Hello Joachim,
How can I acces the files posted in this forums? I always get error 404 Page not found when I try to acces files.

Regards Michael

christopher
Offline
Administrator

2 years 2 months ago #4508 by christopher

Replied by christopher on topic Force between magnet and coil (virtual work)

Hi, I think the link should be fixed now.

Time to create page: 0.122 seconds