# Force between magnet and coil (virtual work)

4 years 6 months ago - 4 years 6 months ago #1173 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: 4 years 6 months ago by Domi.

4 years 6 months ago #1174 by joachim
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:
The following user(s) said Thank You: creativeworker

4 years 5 months ago #1210 by Domi
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:

5 months 3 weeks ago #4505 by Michi
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