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