Download this testcase.

# MigFlow - Copyright (C) <2010-2026>
# <Universite catholique de Louvain (UCL), Belgium
#  Universite de Montpellier, France>
#
# List of the contributors to the development of MigFlow: see AUTHORS file.
# Description and complete License: see LICENSE file.
#
# This program (MigFlow) is free software:
# you can redistribute it and/or modify it under the terms of the GNU Lesser General
# Public License as published by the Free Software Foundation, either version
# 3 of the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program (see COPYING and COPYING.LESSER files).  If not,
# see <http://www.gnu.org/licenses/>.

#!/usr/bin/env python

2D Particle Transport in Vortex Flows

This example simulates the transport of small granular particles by a two-dimensional vortex flow generated around fixed cylindrical obstacles. The coupled fluid–grain problem is solved using MigFlow’s FEM–DEM framework.

Keywords

FEM, DEM, Vortex, Particle Transport

import gmsh
from migflow import fluid, scontact, time_integration
import numpy as np
import os, sys, shutil, random

Output Directory

Create a clean output directory for simulation results.

outputdir = "output" if len(sys.argv) < 2 else sys.argv[1]
shutil.rmtree(outputdir, ignore_errors=True)
os.makedirs(outputdir)

Physical parameters

g = np.array([0, -9.81])  # gravity
rho = 1000  # fluid density
nu = 5e-5  # kinematic viscosity
rhop = 2500  # grains density
r = 2e-5  # grains mean radius

Geometrical parameters

r_in = 2.5e-4
# Rayon des cercles
lc_in = r_in / 2
# Taille de maille aux frontieres interieures
lc_out = r_in
# Taille de maille aux frontieres exterieures
H = 2.6e-3 + 0.005
# Hauteur totale du domaine fluide
L = 8e-3
# Largeur totale du domaine fluide
c = 1.5e-3
# Demi-distance entre les cercles
h = 2e-3
# Hauteur des centres des cercles par rapport a l'origine

Mesh generation

gmsh.initialize()
gmsh.model.add("vortices")
box_tag = gmsh.model.occ.add_rectangle(-L / 2, -H / 2, 0, L, H)
disc0_tag = gmsh.model.occ.add_disk(-c, -H / 2 + h, 0, r_in, r_in)
disc1_tag = gmsh.model.occ.add_disk(+c, -H / 2 + h, 0, r_in, r_in)
gmsh.model.occ.cut([(2, box_tag)], [(2, disc0_tag), (2, disc1_tag)])
gmsh.model.occ.synchronize()


def get_line(x0, x1, y0, y1, eps=1e-6):
    r = gmsh.model.get_entities_in_bounding_box(
        x0 - eps, y0 - eps, -eps, x1 + eps, y1 + eps, eps, 1
    )
    return [tag for dim, tag in r]


gmsh.model.add_physical_group(2, [1], name="domain")
gmsh.model.add_physical_group(1, get_line(-L / 2, +L / 2, +H / 2, +H / 2), name="Top")
gmsh.model.add_physical_group(
    1, get_line(-L / 2, +L / 2, -H / 2, -H / 2), name="Bottom"
)
gmsh.model.add_physical_group(1, get_line(-L / 2, -L / 2, -H / 2, +H / 2), name="Left")
gmsh.model.add_physical_group(1, get_line(+L / 2, +L / 2, -H / 2, +H / 2), name="Right")
gmsh.model.add_physical_group(
    1,
    get_line(-c - r_in, -c + r_in, -H / 2 + h - r_in, -H / 2 + h + r_in),
    name="Disc0",
)
gmsh.model.add_physical_group(
    1,
    get_line(+c - r_in, +c + r_in, -H / 2 + h - r_in, -H / 2 + h + r_in),
    name="Disc1",
)
gmsh.model.mesh.set_size_callback(lambda *args: lc_in)
gmsh.model.mesh.generate(2)

Particle problem initialization

p = scontact.ParticleProblem(2)
bnd_body = p.add_body((0, 0), 0, 0)
p.load_mesh_to_body(bnd_body, None)

h_p = h - r_in
bnd_deposit = np.array(
    [
        (-L / 2, -H / 2 + h_p),
        (-L / 2, -H / 2),
        (+L / 2, -H / 2),
        (+L / 2, -H / 2 + h_p),
    ]
)
n_p_max = 1e4
radii = np.random.uniform(0.7 * r, r, int(n_p_max))
print("Generating deposit particle positions...")
xp, rp = scontact.fastdeposit_2d(bnd_deposit, radii, "deposit.svg")
for xi, ri in zip(xp, rp):
    p.add_particle((xi[0], xi[1]), ri, ri**2 * np.pi * rhop)

Numerical parameters

dt = 5e-3  # time step
tEnd = 3  # final time
outf = 10
t = 0
i = 0


# Function defining the rotation of the inner boundaries
# Parameters:
#              -x:     coordinates of the point on the circle (x for vertical velocity and y for horizontal velocity)
#              -c:     centre of the circle (go along with x)
#              -A:     amplitude of the periodic signal
#              -T:     periode
#              -s:     sign for the rotation
def Bnd(x, c, A, T, s):
    return s * A * (x[:] - c) * np.sin(t * np.pi * 2.0 / T)


yc = -H / 2 + h  # absolute y-centre
v_tang = 2e-2  # max absolute tangential velocity
omega = v_tang / r_in  # max absolute angular velocity
A = omega  # amplitude of periodic signal
T = 1  # period

Fluid problem

fluid = fluid.FluidProblem(2, g, nu * rho, rho)
fluid.load_msh(None)
walls = ["Top", "Bottom", "Left", "Right", "Disc0", "Disc1"]
fluid.set_wall_boundary(walls)
fluid.set_mean_pressure(0)
fluid.set_strong_boundary("Disc1", velocity_x=lambda x: Bnd(x[:, 1], +yc, A, T, +1))
fluid.set_strong_boundary("Disc1", velocity_y=lambda x: Bnd(x[:, 0], +c, A, T, -1))
fluid.set_strong_boundary("Disc0", velocity_x=lambda x: Bnd(x[:, 1], +yc, A, T, -1))
fluid.set_strong_boundary("Disc0", velocity_y=lambda x: Bnd(x[:, 0], -c, A, T, +1))

Computation loop

mass = p.volume() * rhop
while t < tEnd:
    print(f"t/tEnd={t:1.4f}/{tEnd:1.4f}, i={i:1d}")

    fluid.set_particles(
        p.delassus(),
        p.volume(),
        p.position(),
        p.velocity(),
        p.omega(),
        p.contact_forces(),
    )
    if i % outf == 0:
        fluid.write_mig(outputdir, t)
        p.write_mig(outputdir, t)

    fluid.implicit_euler(dt)
    forces = mass * g + fluid.compute_node_force()
    time_integration._advance_particles(p, forces, dt, 1, 1e-3 * r)

    t += dt
    i += 1

Plot

python3 -m migflow.plot.migplot output --actors fluid particles --fluid-field velocity