Force between magnet and coil (virtual work)

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!

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()
    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():

    gfu = GridFunction(fes)

    inv = CGSolver(a.mat, c.mat) = 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)

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


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.


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

Hi, I think the link should be fixed now.