Flow past moving circle cannot implement ALE.move

1

Hi all,

Taking inspiration from the flow past a cylinder tutorial, I’m trying to allow another cylinder initialised at a different place in the domain to move through the fluid. I saw Jorgen’s answer for how to accomplish with constant Dirichlet boundary conditions, although I’m a bit stuck on how to implement this with the inflow profile as specified in the ‘Flow past a cylinder tutorial’. This is the code so far:

# Create mesh
channel = Rectangle(Point(0, 0), Point(2.2, 0.41))
cylinder = Circle(Point(0.2, 0.2), 0.05)
domain = channel - cylinder - following_fish_circle
original_mesh = generate_mesh(domain, 64)

# Define function spaces
V = VectorFunctionSpace(original_mesh, 'CG', 2)
Q = FunctionSpace(original_mesh, 'CG', 1)
# Define boundary and set the boundaries number.
class inflow(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0) and on_boundary

class outflow(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 2.2) and on_boundary

class walls(SubDomain):
    def inside(self, x, on_boundary):
        return (near(x[1], 0) or near(x[1], 0.41)) and on_boundary

class cylinder(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] >= 0.1 and x[0] <= 0.3 and x[1] >= 0.1 and x[1] <= 0.3 and on_boundary

class following_fish(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] >= 1.1 and x[0] <= 1.6 and x[1] >= 0.2 and x[1] <= 0.4 and on_boundary

boundaries = MeshFunction('size_t', original_mesh, 1, original_mesh.domains())
boundaries.set_all(0)
inflow().mark(boundaries, 1)
walls().mark(boundaries, 2)
cylinder().mark(boundaries, 3)
following_fish().mark(boundaries, 4)
outflow().mark(boundaries, 5)

# Define inflow profile
inflow_profile = ('4.0*1.5*x[1]*(0.41 - x[1]) / pow(0.41, 2)', '0')
bcu_inflow = DirichletBC(V, Expression(inflow_profile, degree=2), boundaries, 1)

def boundary_condition_changes(displacement, mesh, bcu_inflow):
    # Define boundary conditions
    bcu_walls = DirichletBC(V, Constant((0, 0)), boundaries, 2)
    bcu_cylinder = DirichletBC(V, Constant((0, 0)), boundaries, 3)
    bcu_following_fish = DirichletBC(V, Constant((-displacement, 0)), boundaries, 4)
    bcp_outflow = DirichletBC(Q, Constant(0), boundaries, 5)
    bcu = [bcu_inflow, bcu_walls, bcu_cylinder, bcu_following_fish]
    bcp = [bcp_outflow]

    if displacement != 0:
        #disp_bc1 = bcu[0]
        disp_bc2 = bcu[1]
        disp_bc3 = bcu[2]
        disp_bc4 = bcu[3]
        #disp_bc5 = bcp[0]
        b_disp = Function(V)
        #disp_bc1.apply(b_disp.vector())
        disp_bc2.apply(b_disp.vector())
        disp_bc3.apply(b_disp.vector())
        disp_bc4.apply(b_disp.vector())
        #disp_bc5.apply(b_disp.vector())

        # Define new mesh
        new_mesh = Mesh(mesh)
        new_boundaries = MeshFunction("size_t", new_mesh, 1)
        new_boundaries.set_values(boundaries.array())
        ALE.move(mesh, b_disp)

        for j in range(700):
            mesh.smooth()

    return bcu, bcp, mesh

This function works without an error when called later in the for loop

for n in range(num_steps):
    if n!=0:
        bcu, bcp, changed_mesh = boundary_condition_changes(displacement=0.01, mesh=changed_mesh, bcu_inflow=bcu_inflow)
        L1, L2, L3, a1, a2, a3, u_, p_, u_n, p_n = initialise_solver(MU=mu, RHO=rho, MESH=changed_mesh)
        A1, A2, A3 = flow_solver(BCU=bcu, BCP=bcp, a1=a1, a2=a2, a3=a3)

Here, initialise_solver and flow_solver are inspired from the ‘flow behind the cylinder’ tutorial. At each time step, the circle at the far end, gets displaced by 0.01 in the negative x direction. The code works but outputs the fluid going through the circle as shown below, instead of around the circle boundary.

velocity_5

I’m a beginner at Fenics and DOLFIN. Could I have some help please?

Thanks!

2

I would start by looking at this mesh function in for instance Paraview;

File("BoundaryMarkers.pvd")<<boundaries

To make sure that you have tagged the cylinder with the appropriate boundary conditions.

Secondly, If you displace your mesh and use smoothing, eventually the mesh will degenerate.
Then remeshing is usually required.
See: https://www.sciencedirect.com/science/article/pii/S0045782520303145
for another approach where one uses a hierarchy of meshes that can move independelty of each other.

3

Hi Jorgen,

It was a visualisation issue, and not a CFD issue. Your tip of looking if the circle had been tagged with the appropriate boundary conditions was spot on. I just applied

plt.clf()

after every u_ and p_ plot and it works very well.

Thanks :slight_smile:

4

Hi Jorgen,

With larger displacements, the mesh does begin to degenerate. Is it possible to apply a simple remeshing with my example. I’m unfamiliar how to go about it?

Instead of the original_mesh being passed here, should I create a new mesh after a threshold cumulative displacement value? How can I ensure that the velocity and pressure vectors computed at the previous time step are carried over when a completely new mesh is generated.

Thanks so much for your time!

Regards,
Vedang

5

Remeshing is not straightforward, as you already mentioned. Let’s say you are able to remesh your geometry after a certain number of time steps. Then you need to interpolate from one grid to the other.
You can for instance use fenicstools ( GitHub - mikaem/fenicstools); see Interpolation nonmatching mesh in parallel · mikaem/fenicstools Wiki · GitHub
to interpolate from one mesh to the other.

6

I’ve had a look at the tutorial, and I’m a little confused. The tutorial demonstrates how to interpolate a function to a single FunctionSpace. But my velocity and pressure solutions depend on a Vector Function Space and another Function Space. How can this interpolation be achieved as the interpolate_nonmatching_mesh_any function does not take two Function Space inputs?

Am I missing something here? Is this interpolation between solution and a Function Space injective? Do you need to interpolate the velocity solution to both Function spaces individually using the above function and do something similar with the pressure solution?