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

Falling particles in a sandglass

This example illustrates the fall of dense solid particles in a viscous fluid under gravity in a two-dimensional hourglass.

Keywords

DEM, FEM

Description

The test demonstrates: - How to initialize a deposit of circular particles. - How to couple the particle solver (scontact) and the fluid solver (fluid).

from migflow import fluid, scontact, time_integration
import numpy as np
import os, shutil, subprocess, sys
import 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)

Geometrical parameters and mesh generation

y0_part = 40e-3  # vertical center of the grains
ly_part = 20e-3  # semi-height of the grains
lx_part = 0.015  # semi-width of the grains
meshfile = "2d_mesh"
subprocess.run(["gmsh", "-2", meshfile + ".geo", "-o", meshfile + ".msh"])

Physical Parameters

g = np.array([0, -9.81])  # gravity
rho = 1  # fluid density
rhop = 1500  # grains density
nu = 1e-5  # kinematic viscosity
r = 1e-4  # grains radius

Simulation Parameters

dt = 2e-4  # time step
tEnd = 1  # final time
outf = 10  # number of iterations between output files
alpha = 0.8  # 2d-3d correction factor for the fluid-particle volume coupling
use_lmgc90 = False  # enable the use of lmgc90

Particle Problem Initialization

p = scontact.ParticleProblem(2)

# Loading of the mesh.msh file specifying physical boundaries name
p.load_msh_boundaries(f"{meshfile}.msh", ["Top", "Box"], material="Material1")

# Generation of the particles in a rectangular grid
x = np.arange(r, lx_part, 2 * r)
y = np.arange(y0_part + r, y0_part + ly_part, 2 * r)
for i in range(x.shape[0]):
    for j in range(y.shape[0]):
        rr = r * (0.95 + random.random() * 0.05)
        p.add_particle((x[i], y[j]), rr, rr**2 * np.pi * rhop, "Material2")

if use_lmgc90:
    from migflow import lmgc90Interface
    print("Deprecated: lmgc90 interface is no longer supported.")
    pass
    # friction=0.1                                     # friction coefficient
    # lmgc90Interface.scontactTolmgc90(outputdir, 2, 0, friction)
    # p = lmgc90Interface.ParticleProblem(2)
else:
    p.set_friction_coefficient(0.4, "Material2", "Material2")  # between particles
    p.set_friction_coefficient(
        0.5, "Material1", "Material2"
    )  # between particles and boundaries

Fluid Problem Initialization

f = fluid.FluidProblem(2, g, nu * rho, rho)
f.load_msh(f"{meshfile}.msh")
f.set_wall_boundary("Top", velocity=[0, 0])
f.set_wall_boundary("Box", velocity=[0, 0])
f.set_mean_pressure(0)

Simulation Loop

Time integration of coupled fluid–particle motion.

t = 0
i = 0
while t < tEnd:
    print(f"{i:4d}, {t:.6g}/{tEnd:.6g}")
    f.set_particles(
        p.delassus(),
        alpha * p.volume(),
        p.position() + p.velocity() * dt,
        p.velocity(),
        p.omega() * 0,
        p.contact_forces(),
    )
    if i % outf == 0:
        p.write_mig(outputdir, t)
        f.write_mig(outputdir, t)
    f.implicit_euler(dt)
    forces = rhop * p.r() ** 2 * np.pi * g + f.compute_node_force()
    time_integration._advance_particles(p, forces, dt, min_nsub=1, contact_tol=1e-3 * r)
    t += dt
    i += 1

Plot

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