Simulate permanent magnet(s) with periodic symmetry along one axis

General
1

Please forgive my inexperience in this field. I’m trying to setup a workflow to simulate magnet arrays (brushless motor rotors) to predict how (small) geometric changes will affect the field shape and strength. I’d like to make use of periodic (rotational) symmetry to decrease setup time / effort.

I have a simple test setup based on one of the site examples. This setup generates fields lines which look good, but as soon as I cut the model in a pie-shape and tag the cut faces as symmetry face pairs the simulation results are inaccurate:

symm
symm1282×564 118 KB

The bit of code used to cut-up the model and ‘tag’ the appropriate faces:

def slice_pie(shape, axis = occ.Z, count = 6, start = 30.0):
    bb = shape.bounding_box
    
    pmin = bb[0]
    pmax = bb[1]

    direction = occ.Pnt(axis.z, axis.x, axis.y)

    center = occ.Pnt(\
            (pmin.x + pmax.x) / 2, \
            (pmin.y + pmax.y) / 2, \
            (pmin.z + pmax.z) / 2 \
        )

    size = max(pmax.x - pmin.x, pmax.y - pmin.y, pmax.z - pmin.z)

    disp = occ.Vec( \
            direction.x * size - center.x, \
            direction.y * size - center.y, \
            direction.z * size - center.z  \
        )
    
    rot_axis = occ.Axis((0,0,0), axis)

    bb1 = bb2 = occ.Box(\
            ( -size, -size, -size ),
            ( size, size, size )
        ).Move(disp).Rotate(rot_axis, 180 + start)
    
    bb2 = bb2.Rotate(rot_axis, 180 - 360 / count)
    
    sliced = shape - bb1
    
    i = 0
    
    pairs = []

    for face in sliced.faces:
        if face.name is None:
            pairs.append(face)
            face.name = f"A_{i}"
            i += 1

    
    sliced = sliced - bb2

    i = 0

    for face in sliced.faces:
        if face.name is None:
            pairs[i] = (pairs[i], face)
            face.name = f"B_{i}"
            i += 1

    return (sliced, pairs)

shape, pairs = slice_pie(shape, magnet_axis, magnet_sections)

for pair in pairs:
    pair[1].Identify(pair[0], "periodic", occ.IdentificationType.PERIODIC) # tried this in reverse order too

geo = occ.OCCGeometry(shape)

The code for calculating the magnetic field:

magnet_mask = ng.CoefficientFunction([0 if mat == 'magnet' else 1 for mat in mesh.GetMaterials()])

fes = ng.HCurl(mesh, order=3, dirichlet="outer", noGrads=True)
print ("ndof =", fes.ndof)
u = fes.TrialFunction()
v = fes.TestFunction()

mu0 = 4 * pi * 1e-7
mur = mesh.MaterialCF({"magnet" : 1000}, default=1) # What unit is used here? Gauss?
nu = 1 / (mu0 * mur)

a = ng.BilinearForm(fes)
a += nu * ng.curl(u) * ng.curl(v) * ng.dx + 1e-8 * nu * u * v * ng.dx

c = ng.Preconditioner(a, "bddc")

f = ng.LinearForm(fes)
mag = mesh.MaterialCF({"magnet" : (magnet_axis.x, magnet_axis.y, magnet_axis.z)}, default=(0,0,0)) # What exactly is (0,0,1) here? The normal for the magnetization direction?
f += mag * ng.curl(v) * ng.dx("magnet")

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

gfu = ng.GridFunction(fes)

with ng.TaskManager():
    solution = ng.solvers.CG(sol=gfu.vec, rhs=f.vec, mat=a.mat, pre=c.mat, maxsteps=256, printrates=True)

ng.Draw (ng.curl(gfu) * magnet_mask, mesh, "B-field")

I suppose I’m missing something here to account for the periodic mesh. How can I make this work as intended?

2

Hi,
first for a rotational identification you need to give the transformation between the faces in the Identify call:

trafo = occ.gp_Trsf.Rotation(Axis(center, axis_vector), angle)
fist_side_faces.Identify(second_side_faces, "periodic", occ.IdentificationType.PERIODIC, trafo)

with this you can also identify all faces in one call (netgen will find the matching faces using the trafo).

Then check in the Netgen gui with View->Viewing Options->Mesh->Show identified points if your identification works. There should be thick violet lines between the identified points.

Last you need to wrap your finite element space with the periodic wrapper:

fes = ng.Periodic(ng.HCurl(....))

When testing this I found a bug in Netgen in the identification for these cake-geometries. They are fixed in the upcoming nightly.

3

Thank you for the hints. I’ve built netgen and ngsolve fresh from git master.

The updated part from the slicing function looks like this now:

    lista = []
    listb = []

    sliced = shape - bb1    
    
    i = 0

    for face in sliced.faces:
        if face.name is None:
            lista.append(face)
            face.name = f"A_{i}"
            i += 1

    
    sliced = sliced - bb2    
    b_faces = sliced.faces;

    i = 0

    for face in sliced.faces:
        if face.name is None:
            listb.append(face)
            face.name = f"B_{i}"
            i += 1

    return (sliced, occ.ListOfShapes(lista), occ.ListOfShapes(listb))

and identifying the faces like this:

    # symmetry = 6
    shape, first_side_faces, second_side_faces = slice_pie(shape, simulation_origin, simulation_axis, symmetry, 0)
    
    trafo = occ.gp_Trsf.Rotation(occ.Axis(simulation_origin, simulation_axis), 360 / symmetry)

    first_side_faces.Identify(second_side_faces, "periodic", occ.IdentificationType.PERIODIC, trafo)

I’ve tried swapping the face sets, using radians as rotation, setting the rotation negative and pretty much any permutation of that, but no violet lines showing things work as intended (results are the same as before too).

What am I (still) doing wrong here?

4

without running the code it’s a bit hard to tell. can you give a minimal example? also possible per dm/email

5

I’m too new a user to upload the file, but the whole script is pretty minimal still:

from math import pi
import ngsolve as ng
from netgen import occ


def add_air_cyl(shape, axis = occ.Z, scale = 5.0):
    bb = shape.bounding_box
    
    pmin = bb[0]
    pmax = bb[1]
    
    center = occ.Pnt(\
            (pmin.x + pmax.x) / 2, \
            (pmin.y + pmax.y) / 2, \
            (pmin.z + pmax.z) / 2 \
        )

    size = max(pmax.x - pmin.x, pmax.y - pmin.y, pmax.z - pmin.z) * scale
    
    cyl = occ.Cylinder((center.x - axis.x * size/2, center.y - axis.y * size/2, center.z - axis.z * size/2), axis, r=size/2, h=size)
    
    cyl.faces.name = "outer"
    
    air = cyl - shape
    air.mat("air")
    
    with_air = occ.Glue([air,shape])

    return with_air
    

def slice_pie(shape, origin = occ.Pnt(0,0,0), axis = occ.Z, count = 6, start = 30.0):
    bb = shape.bounding_box
    
    pmin = bb[0]
    pmax = bb[1]

    direction = occ.Pnt(axis.z, axis.x, axis.y)

    center = occ.Pnt(\
            (pmin.x + pmax.x) / 2, \
            (pmin.y + pmax.y) / 2, \
            (pmin.z + pmax.z) / 2 \
        )

    size = max(pmax.x - pmin.x, pmax.y - pmin.y, pmax.z - pmin.z)

    disp = occ.Vec( \
            direction.x * size - center.x, \
            direction.y * size - center.y, \
            direction.z * size - center.z  \
        )
    
    rot_axis = occ.Axis(origin, axis)

    bb1 = bb2 = occ.Box(\
            ( -size, -size, -size ),
            ( size, size, size )
        ).Move(disp).Rotate(rot_axis, 180 + start)
    
    bb2 = bb2.Rotate(rot_axis, 180 - 360 / count)
    
    lista = []
    listb = []

    sliced = shape - bb1    
    
    i = 0

    for face in sliced.faces:
        if face.name is None:
            lista.append(face)
            face.name = f"A_{i}"
            i += 1

    
    sliced = sliced - bb2    
    b_faces = sliced.faces;

    i = 0

    for face in sliced.faces:
        if face.name is None:
            listb.append(face)
            face.name = f"B_{i}"
            i += 1

    return (sliced, occ.ListOfShapes(lista), occ.ListOfShapes(listb))
    

def make_magnet(axis = occ.Z, r = 0.3, h = 2.0, resolution = 0.1):
    if not hasattr(make_magnet, "counter"):
        make_magnet.counter = 0
    matname = f"magnet_{make_magnet.counter}"
    make_magnet.counter += 1

    shape = occ.Cylinder((-0.5 * h * axis.x, -0.5 * h * axis.y, -0.5 * h * axis.z), axis, r, h)
    shape.mat(matname)
    shape.maxh = resolution
    shape.faces.name = "magnet"
    shape.faces.col = (1,0,0)
    
    yield shape;
    
    global mesh
    mur = mesh.MaterialCF({"magnet" : 1000}, default=1) # What unit is used here? Gauss?
    nu = 1 / (mu0 * mur)
    yield nu

    mag = mesh.MaterialCF({matname : (axis.x, axis.y, axis.z)}, default=(0,0,0))
    yield mag * ng.curl(v) * ng.dx(matname)

ng.ngsglobals.msg_level = 3

mu0 = 4 * pi * 1e-7

resolution = 2
symmetry = 6
simulation_axis = occ.Z
simulation_origin = occ.Pnt(0,0,0)
magnet_axis = occ.Z 
magnet_diameter = 6.0 #mm
magnet_length = 10.0 #mm
air_scale = 5.0

magnet_resolution = min(magnet_diameter, magnet_length) / (15 * resolution)
overall_resolution = air_scale * min(magnet_diameter, magnet_length) / (5 * resolution)

print ("Resolution: ", overall_resolution, "/", magnet_resolution)

mag_gen = make_magnet(magnet_axis, r=magnet_diameter/2, h=magnet_length, resolution=magnet_resolution)

magnet = next(mag_gen)

shape = add_air_cyl(magnet, magnet_axis)

if symmetry > 1:
    shape, first_side_faces, second_side_faces = slice_pie(shape, simulation_origin, simulation_axis, symmetry, 0)
    
    # OpenCascade takes radians (from what I've seen in the source), but ngsolve seems to use degrees and do a translation?
    # I tried 2*pi instead of 360, also both negative
    trafo = occ.gp_Trsf.Rotation(occ.Axis(simulation_origin, simulation_axis), 360 / symmetry)

    first_side_faces.Identify(second_side_faces, "periodic", occ.IdentificationType.PERIODIC, trafo)

    # for pair in pairs:
    #     pair[1].Identify(pair[0], "periodic", occ.IdentificationType.PERIODIC, trafo)
    #     print("pair: ", pair[0].name, " - ", pair[1].name)

geo = occ.OCCGeometry(shape)

ng.Draw (geo)

nmesh = geo.GenerateMesh(maxh=overall_resolution, curvaturesafety=1)

mesh = ng.Mesh(nmesh).Curve(3)

mesh.GetMaterials()
mesh.GetBoundaries()

ng.Draw (mesh)

magnet_mask = ng.CoefficientFunction([0 if mat == 'magnet' else 1 for mat in mesh.GetMaterials()])

fes = ng.HCurl(mesh, order=3, dirichlet="outer", noGrads=True)
if symmetry > 0:
    fes = ng.Periodic(fes)
    
print ("ndof =", fes.ndof)
u = fes.TrialFunction()
v = fes.TestFunction()

a = ng.BilinearForm(fes)
nu = next(mag_gen)
a += nu * ng.curl(u) * ng.curl(v) * ng.dx + 1e-8 * nu * u * v * ng.dx

c = ng.Preconditioner(a, "bddc")

f = ng.LinearForm(fes)
f += next(mag_gen)

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

gfu = ng.GridFunction(fes)

with ng.TaskManager():
    solution = ng.solvers.CG(sol=gfu.vec, rhs=f.vec, mat=a.mat, pre=c.mat, maxsteps=256, printrates=True)

ng.Draw (ng.curl(gfu) * magnet_mask, mesh, "B-field")
# ng.Draw (ng.grad(gfu) * magnet_mask, mesh, "B-field")
# ng.Draw (filtered_b_field, mesh, "B-field")
# ng.Draw (b_field, mesh, "B-field")

# vtk = ng.VTKOutput(mesh, coefs=[ng.curl(gfu)], names=["B-field"], filename="bfield", subdivision=2)

# vtk.Do()
# vtk.Do(vb=ng.BND)
6

The problem is that after cutting the second time the faces that you have kept in lista are not in the geometry any more (they are also cut again). That’s why we like to work with string querying of the face names there, this works:

shape.faces["B_.*"].Identify(shape.faces["A_.*"], "periodic", occ.IdentificationType.PERIODIC, trafo)

upside is also you do not have to keep list around. Just name everything correctly and use the names afterwards. You can check the names in the gui by double clicking on the mesh elements.

You can use regular expressions like "B_.*" to select all faces that start with “B_”.

7

Oh that’s much easier indeed, thank you!

I checked and with your change it works beautifully. Thank you again for all the help!

8

Last follow-up question on this, if you don’t mind: Is there a way to visualize the mesh and results as if not using symmetry (without solving it again)?

I imagine some way to revolve the mesh and results around the simulation axis n-times? I think I can get the mesh done, but I’m not sure how to do the same for the results.

9

Hm in ngsolve not really, although with the new webgui I think it might be possible with a bit of hacking.
But as far as I know paraview supports this so you could export the vtu and do it there.

10

Thanks! I’ll give paraview a try for that.

For those playing with above code, I made one (imho) worthwhile addition to the slice_pie function, just before the return statement add:

    for subshape in sliced.solids:
        print ("subshape maxh:", subshape.maxh)
        for face in subshape.faces["[AB]_[0-9]+"]:
            face.maxh = min(face.maxh, subshape.maxh)

This makes the slice faces inherit the resolution of the sliced solids, for more consistent meshes.