Sector periodicity for rotating domains

dolfinx
1

I am considering a rotating electric machine which is efficiently solved for periodic sectors specially in 3D. Below is an exemplary image to convey the idea:

image
image1024×768 14.2 KB

I was wondering whether we have tools in dolfinx that could be used for enforcing periodicity and continuity required by these problems? I have used dolfinx_mpc for enforcing continuity across non-conformal interfaces but it throws error as soon as the interfaces slide through each other. I would appreciate any advice (even negative ones) from informed users and developers.

Thanks

2

There is a flag in DOLFINx_MPC that can be used to allow partial contact/periodicity: (allow_missing_masters=True) as seen in: dolfinx_mpc/python/src/dolfinx_mpc/multipointconstraint.py at 0894ca6368640df9ad53e4b7aa57737208d6ea4f · jorgensd/dolfinx_mpc · GitHub

Of course you would have to recompute the periodicity as the mesh moves and the periodic conditions changes.

3

Thanks for the pointer! As I look at it, it appears to be a remedy for two surfaces partially out of contact due to movement. The segment still in contact will connect while the one out of touch will be subject to natural boundary condition (right?). But in my case, the out of contact boundary is still in contact with the neighboring sector not explicitly modeled due to sector periodicity. In fact, even if it has no visible contact, it must still have full contact. Based on the docstring, I see no possibility to specify this modulo periodic aspect of the problem. Could you confirm my understanding?

4

Wouldn’t that be posing a sequence of constraints on the interface, where each of them allow_missing_masters?

5

I am not that familiar to share thoughts on your suggestion. If I provide you an MWE, could you use it to try it out?

#%% input parameters
fov_R = 1 # radius of the field of view
R_int = 0.5 # radial position of the interface
n_stat_poles = 8 # number of poles i.e. sector
sec_angle = 2*np.pi/n_stat_poles
rot_angle = 350 * np.pi / 180 # angle of rotation of the rotor
fem_order = 2
rot_dom_idx = (1, ) # rotating (top) domain, by inspection
stat_dom_idx = (2, 3)
mag_dom = (2,) # by inspection
#%% cad model

gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)
gmsh.clear()

occ = gmsh.model.occ
gdim = 2

sec_angle = 2*np.pi/n_stat_poles

c1 = occ.addCircle(0, 0, 0, fov_R, angle1=np.pi/2-sec_angle/2, angle2=np.pi/2+sec_angle/2)
pnts = occ.getEntities(0)
p0 = occ.addPoint(0, 0, 0)
l1 = occ.addLine(p0, pnts[0][1])
l2 = occ.addLine(p0, pnts[1][1])
loop = occ.addCurveLoop([c1, l1, l2])
slice_ = occ.addPlaneSurface([loop])
slice_copy = occ.copy([(gdim, slice_)])

d1 = occ.addDisk(0, 0, 0, fov_R, fov_R)
d2 = occ.addDisk(0, 0, 0, R_int, R_int)
rot = occ.cut([(gdim, d1)], [(gdim, d2)])
stat = occ.addDisk(0, 0, 0, R_int, R_int)

occ.synchronize()
rot = occ.intersect([(gdim, slice_)], rot[0])

stat = occ.intersect(slice_copy, [(gdim, stat)])
r1 = occ.addRectangle(-0.025, 0.25, 0, 0.05, 0.2)
stat = occ.fragment(stat[0], [(gdim, r1)])

occ.synchronize()

all_doms = occ.getEntities(gdim)

all_doms = gmsh.model.getEntities(gdim)
for j, dom in enumerate(all_doms):
    gmsh.model.addPhysicalGroup(dom[0], [dom[1]], j + 1)  # create the main group/node

# number all boundaries
all_edges = gmsh.model.getEntities(gdim - 1)
for j, edge in enumerate(all_edges):
    gmsh.model.addPhysicalGroup(edge[0], [edge[1]], edge[1])  # create the main group/node

# number all boundaries
all_edges = gmsh.model.getEntities(gdim - 2)
for j, edge in enumerate(all_edges):
    gmsh.model.addPhysicalGroup(edge[0], [edge[1]], edge[1])  # create the main group/node

gmsh.model.mesh.generate(gdim)
gmsh.model.mesh.refine()

mpi_rank = 0
dolfinx_model = dolfinx.io.gmsh.model_to_mesh(gmsh.model, MPI.COMM_WORLD, mpi_rank, gdim)
mesh, ct, ft, rt = dolfinx_model.mesh, dolfinx_model.cell_tags, dolfinx_model.facet_tags, dolfinx_model.ridge_tags

# gmsh.fltk.run()

#%% physical model

def rotate_mesh(mesh, ct, dom_idx, angle):
    '''Rotate cells in submodains identified in dom_idx tuple by the specified angle'''
    # https://fenicsproject.discourse.group/t/how-to-rotate-a-mesh-in-fenics/7162/6

    from scipy.spatial.transform import Rotation
    gdim = 2
    
    # identify vertices that will be rotated
    rot_cell_indices = [ct.find(idx) for idx in dom_idx]
    rot_cell_indices = np.hstack(rot_cell_indices)
    # vertices are shared between cells. Therefore, we keep only unique
    cell_vertex_indices = dolfinx.cpp.mesh.entities_to_geometry(mesh._cpp_object, gdim, rot_cell_indices, False)
    cell_vertex_indices = np.unique(cell_vertex_indices.ravel())
    # finally apply rotation on identified vertices
    rotation = Rotation.from_euler("z", angle, degrees=False)
    mesh.geometry.x[cell_vertex_indices,:] = rotation.apply(mesh.geometry.x[cell_vertex_indices])

# rotate rotor
rotate_mesh(mesh, ct, np.hstack(rot_dom_idx), rot_angle)

Omega_A = dolfinx.fem.functionspace(mesh, ("CG", fem_order))
A, At = ufl.TrialFunction(Omega_A), ufl.TestFunction(Omega_A)
A0 = dolfinx.fem.Function(Omega_A) # for previous time step

A_3d, At_3d = ufl.as_vector((0, 0, A)), ufl.as_vector((0, 0, At))
A0_3d = ufl.as_vector((0, 0, A0))

def curl_2d(A):
    return ufl.as_vector((A[2].dx(1), -A[2].dx(0), 0))

dx = ufl.Measure("dx", domain=mesh, subdomain_data=ct, metadata={"quadrature_degree": 8})
ds = ufl.Measure("ds", domain=mesh, subdomain_data=ft, metadata={"quadrature_degree": 8})
dS = ufl.Measure("dS", domain=mesh, subdomain_data=ft, metadata={"quadrature_degree": 8})

a = ufl.inner(1/1000*curl_2d(A_3d), curl_2d(At_3d))*dx(rot_dom_idx[0])
a += ufl.inner(curl_2d(A_3d), curl_2d(At_3d))*dx(tuple(stat_dom_idx))
# add permanent magnet
L = ufl.inner(dolfinx.fem.Constant(mesh, (0.0, 1.0, 0.0)), curl_2d(At_3d))*dx(mag_dom)

dbc = []

#%% MPC
#  
# https://github.com/jorgensd/dolfinx_mpc/blob/main/python/demos/demo_periodic_geometrical.py

# periodicity constraint
def periodic_boundary_stator(x):
    r = np.sqrt(x[0]**2+x[1]**2)
    theta = np.arctan2(x[1], x[0])
    sector_angle = 2*np.pi/n_stat_poles
    bnd_marker = np.isclose(theta, np.pi/2 - sector_angle/2, atol=tol)
    stator_marker = r <= R_int
    return bnd_marker & stator_marker

def periodic_boundary_rotor(x, rot_angle):
    r = np.sqrt(x[0]**2+x[1]**2)
    theta = np.arctan2(x[1], x[0])
    sector_angle = 2*np.pi/n_stat_poles
    bnd_marker = np.isclose(theta, np.pi/2 - sector_angle/2 + rot_angle, atol=tol)
    rotor_marker = r >= R_int
    return bnd_marker & rotor_marker

#Cyclic Symmetry Condition
def periodic_relation(x, rot_angle=0.0):
    out_x = np.zeros(x.shape)
    sector_angle = 2*np.pi/n_stat_poles
    rot_mat = np.array([
        [np.cos(sector_angle+rot_angle), -np.sin(sector_angle+rot_angle)],
        [np.sin(sector_angle+rot_angle), np.cos(sector_angle+rot_angle)]
    ])
    out_x[:2,:] = rot_mat @ x[:2,:]
    return out_x

mpc = dolfinx_mpc.MultiPointConstraint(Omega_A)
tol = 250 * np.finfo(dolfinx.default_scalar_type).resolution
interface_bnd_idx = (12, 17) # rotor boundary, stator boundary - by visible inspection
mpc.create_contact_inelastic_condition(ft, interface_bnd_idx[0], interface_bnd_idx[1], eps2=tol, allow_missing_masters=True)
# stator periodicity
mpc.create_periodic_constraint_geometrical(Omega_A, periodic_boundary_stator, periodic_relation, dbc, scale=-1.0)
# rotor periodicity
mpc.create_periodic_constraint_geometrical(Omega_A, 
                                           lambda x: periodic_boundary_rotor(x, rot_angle), 
                                           lambda x: periodic_relation(x, rot_angle), dbc, scale=-1.0)
mpc.finalize()

petsc_options = {"pc_type": "lu", "pc_factor_mat_solver_type": "mumps"}
problem = dolfinx_mpc.problem.LinearProblem(a, L, mpc, bcs=dbc, petsc_options=petsc_options)
Az = problem.solve() # Wb/m
Az.x.scatter_forward()

Omega_B = dolfinx.fem.functionspace(mesh, ("DG", fem_order-1, (3,)))
B_stat = dolfinx.fem.Function(Omega_B)
B_expr = dolfinx.fem.Expression(curl_2d(ufl.as_vector((0, 0, Az))), Omega_B.element.interpolation_points)
B_stat.interpolate(B_expr)

#%% Generate pyvista plot

u_topology, u_cell_types, u_geometry = dolfinx.plot.vtk_mesh(Omega_A)

u_grid = pyvista.UnstructuredGrid(u_topology, u_cell_types, u_geometry)
u_grid.point_data["Az"] = Az.x.array
u_grid.set_active_scalars("Az")

u_plotter = pyvista.Plotter()
u_plotter.add_mesh(u_grid, show_edges=False, cmap="viridis")
u_plotter.view_xy()
if not pyvista.OFF_SCREEN:
    u_plotter.show()

# |B|
u_topology, u_cell_types, u_geometry = dolfinx.plot.vtk_mesh(B_stat.function_space)

u_grid = pyvista.UnstructuredGrid(u_topology, u_cell_types, u_geometry)
u_grid.point_data["B"] = np.sqrt(B_stat.x.array[::3]**2 + B_stat.x.array[1::3]**2)
u_grid.set_active_scalars("B")

u_plotter = pyvista.Plotter()
u_plotter.add_mesh(u_grid, show_edges=False, cmap="viridis")
u_plotter.view_xy()
if not pyvista.OFF_SCREEN:
    u_plotter.show()