Download :download:`this testcase <2d_vortices.py>`.

.. code-block:: python
  
  # 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

.. code-block:: python
  
  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.

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

Physical parameters
-------------------

.. code-block:: python
  
  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
----------------------

.. code-block:: python
  
  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
---------------

.. code-block:: python
  
  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
-------------------------------

.. code-block:: python
  
  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
--------------------

.. code-block:: python
  
  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
-------------

.. code-block:: python
  
  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
----------------

.. code-block:: python
  
  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
----
.. code-block:: shell

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

