Moving mesh causes entangled cells

mesh
1

I am trying to implement fluid-rigid body interaction with a gear-like structure moving in fluids and here I am trying to move the mesh around. I modify the code from https://fenicsproject.org/qa/13470/ale-move-class-and-meshes/

from dolfin import *
from mshr import *
import matplotlib as plt

def solve_linear_elasticity(mesh, boundaries, d):
    c = Constant(d)

    V = VectorFunctionSpace(mesh, "Lagrange", 1)
    u = TrialFunction(V)
    v = TestFunction(V)

    E, nu = 10.0, 5.0
    mu = E/(2.0*(1.0 + nu))
    lmbda = E*nu/((1.0 + nu)*(1.0 -2.0*nu))
    sigma = 2*mu*sym(grad(u)) + lmbda*tr(grad(u))*Identity(2)
    F = inner(sigma, grad(v))*dx 
    a, L = lhs(F), rhs(F)

    bcs = [DirichletBC(V, Constant((0.0, 0.0)), boundaries, 1),
       DirichletBC(V.sub(0), c, boundaries, 2)]

    displacement = Function(V)
    solve(a==L, displacement, bcs)
    return displacement

def update_mesh(mesh, displacement, boundaries):

    new_mesh = Mesh(mesh)
    new_boundaries = MeshFunction("size_t", new_mesh, 1)
    new_boundaries.set_values(boundaries.array())
    ALE.move(new_mesh, displacement)
    return new_mesh, new_boundaries


# Original mesh
# Generate mesh
domain_vertices = [Point(20, 23.125),
               Point(20, 21.25),
               Point(17.2925, 21.5625),
               Point(18.9175, 20.625),
               Point(17.2925, 18.4375),
               Point(18.9175, 19.375),
               Point(20, 16.875),
               Point(20, 18.75),
               Point(22.705, 18.4375),
               Point(21.0825, 19.375),
               Point(22.705, 21.5625), 
               Point(20.995, 20.8275)]                 
inclusion = Polygon(domain_vertices)
mesh = generate_mesh(Rectangle(Point(0, 0), Point(40, 40))-inclusion, 10)

class on_inclusion(SubDomain):
def inside(self, x, on_boundary):
    return (on_boundary and not(near(x[1], 0) or near(x[1], 40) or near(x[0], 0) or near(x[0], 40)))

# Initialize mesh function for boundary domains
boundaries = MeshFunction('size_t', mesh, 1, mesh.domains())
boundaries.set_all(1)
inc = on_inclusion()
inc.mark(boundaries, 2)

displacement = solve_linear_elasticity(mesh, boundaries, 5)
mesh, boundaries = update_mesh(mesh, displacement, boundaries)

plt.figure
plot(mesh)

Before the moving it looks like this (the gear structure is in the middle)
original_mesh
After the moving the meshes becomes this
mesh

It seems the cells are entangled. Is there anyway to prevent this?

Besides, I would like to include the rotation of the object. Then the boundary condition would vary at each point on the gear, and I wonder how should I specify it?

1 Like
2

The issue is not moving the mesh, the issue is how you solve your elasticity problem. You have not properly marked the outer boundaries with 1. Right now you have marked every boundary with 1, thus only the vertices of your problem is moved.
The following code corrects this:

class on_inclusion(SubDomain):
    def inside(self, x, on_boundary):
        return (on_boundary and not(near(x[1], 0) or near(x[1], 40) or near(x[0], 0) or near(x[0], 40)))

class all_boundaries(SubDomain):
    def inside(self, x, on_boundary):
        return (on_boundary)


# Initialize mesh function for boundary domains
boundaries = MeshFunction('size_t', mesh, 1, mesh.domains())
boundaries.set_all(0)
all = all_boundaries()
all.mark(boundaries,1)
inc = on_inclusion()
inc.mark(boundaries, 2)

Sidenote:
I find your Lame parameters strange, as nu typically is in the range (0,0.5).

If you want to do multiple revolutions of your geometry, I would suggest using a cut finite element method, using either levelsets or multiple meshes to describe the rotating object.
See for instance this paper: https://msp.org/camcos/2015/10-2/p01.xhtml

2 Likes