Download this example.

Bidimensional fall of a cloud made of particles in a viscous fluid.

This examples illustrates how to define the fluid problem as well as the particles one. All the relevant mesh information is directly loaded from the GMSH api.

import os
import time
import gmsh
import shutil
import ctypes
import numpy as np
from scipy.spatial import KDTree
from migflow import fluid
from migflow import scontact
from migflow import time_integration

use_callback = True

First, let’s create the output directory. Define output directory

outputdir = "output-drop"
shutil.rmtree(outputdir,ignore_errors=True)
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)

Mesh generation

The mesh geometry is created with the python GMSH api.

height = 2                                      # domain height
width = 0.1                                     # domain width
origin = [-width/2, -height/2]                  # leftmost bottom corner

def gen_geometry(width, height, origin=np.array([0,0])):
    origin = np.asarray(origin)
    gmsh.model.add("box")
    gmsh.model.occ.add_rectangle(origin[0], origin[1], 0, width, height)
    gmsh.model.occ.synchronize()
    def get_line(x0, x1, eps =1e-6):
        r = gmsh.model.get_entities_in_bounding_box(x0[0]-eps, x0[1]-eps, -eps, x1[0]+eps, x1[1]+eps, eps, 1)
        return list(tag for dim, tag in r)
    h, w = height, width
    gmsh.model.add_physical_group(1,get_line(origin+[0,0], origin+[w,0]), name="Bottom")
    gmsh.model.add_physical_group(1,get_line(origin+[0,h], origin+[w,h]), name="Top")
    gmsh.model.add_physical_group(1,get_line(origin+[0,0], origin+[0,h])+get_line(origin+[w,0], origin+[w,h]), name="Lateral")
    gmsh.model.add_physical_group(2,[1], name="domain")

The mesh size can be chosen based on three approaches; based on a computed field, based on a callback or based on a uniform mesh size.


The mesh size is computed based on the gradient variation. The field is then given to gmsh which will generate the mesh.

def gen_mesh_field(f):
    lcmin = 0.0015/4
    lcmax = 0.015/8
    grad = f.fields_gradient()
    grad_u, grad_v, grad_p = grad[:,0,:], grad[:,1,:], grad[:,2,:]
    grad_v = np.linalg.norm(grad_v, axis=1)
    grad_v_min = np.min(grad_v)
    grad_v_max = np.max(grad_v)
    lv = (grad_v_max - grad_v_min)/(grad_v-grad_v_min) * lcmin
    size = np.maximum(lv, lcmax)

    x = f.coordinates()[f.elements()]
    size_view = gmsh.view.add("size")
    data = np.c_[x.swapaxes(1,2).reshape(-1,9), size[f.elements()]]
    size_field = gmsh.model.mesh.field.add("PostView")
    gmsh.model.mesh.field.setNumber(size_field, "ViewTag", size_view)
    gmsh.model.mesh.field.setAsBackgroundMesh(size_field)
    gmsh.view.add_list_data(size_view, "ST", data.shape[0], data.reshape(-1))
    gmsh.model.mesh.clear()
    gmsh.model.mesh.generate(2)
    gmsh.view.remove(size_view)

The mesh size is computed based on a callback function. Here a refinement is done close to the particles position.

def gen_mesh_callback(xp):
    tree = KDTree(xp)
    def size_f(x):
        dist, _ = tree.query(x[:,:2])
        distmin = 0.01
        lcmin = 0.0015
        lcmax = 0.015
        alpha = np.clip((dist[0]-distmin)/0.1, 0, 1)
        size = lcmin*(1-alpha) + lcmax*alpha
        return size
    gmsh.model.mesh.set_size_callback(lambda dim, tag, x, y, z, lc: size_f(np.array([[x,y]])))
    gmsh.model.mesh.clear()
    gmsh.model.mesh.generate(2)
    gmsh.model.mesh.remove_size_callback()

The mesh is generated based on a constant mesh size.

def gen_mesh_uniform(size):
    gmsh.model.mesh.clear()
    gmsh.model.mesh.set_size_callback(lambda dim, tag, x, y, z, lc: size)
    gmsh.model.mesh.generate(2)
    gmsh.model.mesh.remove_size_callback()

Particle Problem

The particle problem is created and its dimension is provided.

p = scontact.ParticleProblem(2)

The particle properties are defined. All the particles are assumed to be spherical.

r    = 25e-6                                       # particles radius
rhop = 2450                                         # particles density
compacity = 0.2                                     # solid volume fraction in the drop
rout = 3.3e-3                                       # drop radius

The particle positions are initialised randomly in a circular domain to obtain a given compacity.

def genInitialPosition(p, r, rout, rhop, compacity, origin) :
    """Set all the particles centre positions and create the particles objects

    Keyword arguments:
    p -- Particle Problem
    r -- max radius of the particles
    rout -- outer radius of the cloud
    rhop -- particles density
    compacity -- initial compacity in the cloud
    """
    # Space between the particles to obtain the expected compacity
    N = compacity*(rout/r)**2
    bodies = np.asarray([], int)
    while p.n_particles() < N:
        xyp = np.random.uniform(-rout, rout, 2)
        if np.hypot(xyp[0], xyp[1]) < rout:
            if p.n_particles() == 0:
                d = 1
            else:
                d = np.min(np.linalg.norm(p.position()-xyp[None,:],axis=1))
            if d > 2.1*r:
                body = p.add_particle(xyp, r, r**2 * np.pi * rhop)
                bodies = np.append(bodies,body)
    # Shift of the particles to the top of the box
    p.body_position()[bodies, :] += origin
    return bodies
bodies = genInitialPosition(p, r, rout, rhop, compacity, np.array([0,1.9*height/4]))

Fluid Problem

The fluid is described by its dimension, the external volume force applied, its dynamic viscosity and its density.

g   = np.array([0,-9.81])                           # gravity
rho = 1030                                          # fluid density
nu  = 1.17/(2*rho)                                  # kinematic viscosity
mu  = rho*nu                                        # dynamic viscosity
f   = fluid.FluidProblem(2,g,[nu*rho],[rho], usolid=True)

The mesh is created with a refinment close to the particles position. The boundaries are loaded either from a .msh file or from the current model loaded with the GMSH api (if None as a filename is given).

gen_geometry(width, height, origin)
if use_callback:
    gen_mesh_callback(p.position()[p.r()[:,0]>0])
else:
    gen_mesh_uniform(0.015)
    f.load_msh(None)
p.load_msh_boundaries(None, ["Bottom","Top","Lateral"])
## %%
# The mesh can be loaded through the GMSH api if None is given or through a mesh filename.
# The boundary conditions are described by their physical name.
# To fully describe the pressure, its mean value over all the nodal values is set to zero.
f.load_msh(None)
f.set_wall_boundary("Bottom")
f.set_wall_boundary("Top")
f.set_wall_boundary("Lateral")
f.set_mean_pressure(0)

FEM-DEM coupling

The presence of particles is given to the fluid through the fluid volume fraction, i.e. the porosity and through a drag parametrization. All the relevant informations needed for the fluid is given by :

f.set_particles(p.delassus(), p.volume(), p.position(), p.velocity(), p.omega(), p.contact_forces())

Time integration

The numerical parameters are defined given and the initial conditions are written in the output directory.

outf = 10                                           # number of iterations between output files
remeshing=20                                         # time step to adapt the mesh
dt   = 5e-3                                         # time step
tEnd = 50                                           # final time
t    =  0                                           # initial time
ii   =  0                                           # initial iteration

Write chosen fields into the output

def get_fields(fluid):
    y = fluid.coordinates_fields()[fluid.pressure_index()][:,1]
    y = y.reshape(-1,1)
    p1_element = fluid.get_p1_element()
    return {"pressure": (fluid.pressure(), p1_element),
            "velocity": (fluid.velocity(), p1_element),
            "porosity": (fluid.porosity(), p1_element),
            "dynamic_pressure":(fluid.pressure()-rho*g[1]*y, p1_element)}

p.write_mig(outputdir, t)
f.write_mig(outputdir, t, get_fields(f))

The computational loop can start. The particles phase will use sub-iterations to keep a stable method. The minimal number of subiterations is given by min_nsub. The external forces given to the particles have to be given at each time step. A predictor-corrector method can also be used by setting the flag to True.

tic = time.process_time()
forces = g*np.pi*p.r()**2*rhop
mass = np.pi*p.r()**2*rhop
while t < tEnd :
    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.process_time() - tic))
    f.set_particles(p.delassus(), p.volume(), p.position(), p.velocity(), p.omega(), p.contact_forces())
    tic = time.process_time()
    f.implicit_euler(dt, check_residual_norm=1)
    print("implicit : ", time.process_time()-tic)
    fext = f.compute_node_force(dt) + g*mass
    tic = time.process_time()
    time_integration._advance_particles(p, fext, dt, 1, 1e-3*r)
    print("NSCD : ", time.process_time()-tic)

    # time_integration.iterate(f, p, dt, min_nsub=2, external_particles_forces=forces, use_predictor_corrector=False)
    if ((ii+1)%remeshing==0):
        number_p = f.n_particles
        position_p = f.particle_position()
        volume_p = f.particle_volume()
    if (ii%remeshing==0 and ii != 0):
        if use_callback:
            gen_mesh_callback(p.position()[p.r()[:,0]>0])
        else:
            gen_mesh_field(f)
        f.adapt_mesh(old_n_particles=number_p, old_particle_position=position_p, old_particle_volume=volume_p)
        f.set_particles(p.delassus(), p.volume(), p.position(), p.velocity(), p.omega(), p.contact_forces())
    t += dt
    # Output files writting
    if ii%outf == 0 :
        p.write_mig(outputdir, t)
        f.write_mig(outputdir, t, get_fields(f))
    ii += 1