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

Study of the drag on a fixed disk in a 2D granular flow
=======================================================

This example illustrates how to compute the drag force exerted on
a fixed circular obstacle immersed in a granular flow.
The particles fall under gravity between two lateral walls and impact
the fixed inner cylinder.

Keywords
--------
DEM, Generation, Post-processing

.. code-block:: python
  
  
  
  import os, sys, shutil, subprocess
  import numpy as np
  from migflow import scontact, time_integration
  

Output Directory
----------------
The output directory is (re)created

.. code-block:: python
  
  outputdir = "output" if len(sys.argv) < 2 else sys.argv[1]
  geomesh_filename = "2d_mesh.geo"
  mesh_filename = f"{outputdir}/mesh.msh"
  csv_file = f"{outputdir}/Drag.csv"
  shutil.rmtree(outputdir, ignore_errors=True)
  os.makedirs(outputdir)
  subprocess.call(["gmsh", "-2", geomesh_filename, "-o", mesh_filename])
  

Geometrical parameters
----------------------

.. code-block:: python
  
  r = 6e-3  # particle radius [m]
  h, w = 1.0, 0.28  # particle region height & width [m]
  h0 = -0.5  # particle bed offset
  r_outer = 0.025  # radius of the fixed disk [m]
  

Physical Parameters
-------------------

.. code-block:: python
  
  rhop = 2500  # particle density [kg/m³]
  g = np.array([0, -9.81])  # gravity [m/s²]
  friction = 0.2  # friction coefficient
  

Particle Problem Initialization
-------------------------------

Particle problem setup
----------------------
Load mesh boundaries (requires a `mesh.msh` file)

.. code-block:: python
  
  p = scontact.ParticleProblem(2)
  p.load_msh_boundaries(mesh_filename, ["Inner"], material="Steel")
  p.load_msh_boundaries(
      mesh_filename,
      ["Right", "Left", "RightUp", "RightDown", "LeftUp", "LeftDown"],
      material="Plexi",
  )
  p.set_friction_coefficient(friction, "Sand", "Sand")
  p.set_friction_coefficient(friction, "Sand", "Steel")
  
  

Particle initialization
-----------------------
Add granular particles around the fixed disk.

.. code-block:: python
  
  def gen_rect(origin, w, h, step):
      """Generate coordinates on a rectangular grid.
  
      Parameters
      ----------
      origin : list(float)
          Origin of the top-left corner.
      w : float
          Width of the rectangular domain.
      h : float
          Height of the rectangular domain.
      step : float
          Particle radius (spacing control).
  
      Returns
      -------
      tuple of np.ndarray
          Arrays of x and y coordinates.
      """
      eps = 1e-8
      x = np.arange(-w / 2 + step - eps, w / 2 - step + eps, 2 * step) + origin[0]
      y = np.arange(step, h - step, 2 * step) + origin[1]
      x, y = np.meshgrid(x, y)
      return x.ravel(), y.ravel()
  
  
  # Fill the domain around the disk
  x, y = gen_rect([0, h0], w, 1.5 * h, r)
  for xi, yi in zip(x, y):
      if xi**2 + yi**2 > (r_outer + r) ** 2:
          p.add_particle((xi, yi), r, np.pi * r**2 * rhop, "Sand")
  # Add a reservoir of particles above
  x, y = gen_rect([0, -2 * h0], 2 * w, 0.5 * h, r)
  for xi, yi in zip(x, y):
      p.add_particle((xi, yi), r, np.pi * r**2 * rhop, "Sand")
  

Numerical Parameters
--------------------

.. code-block:: python
  
  dt = 1e-3  # time step [s]
  tEnd = 2.0  # final time [s]
  outf = 10  # output frequency
  t, i = 0.0, 0
  
  

Helper function to accumulate the drag forces
---------------------------------------------

.. code-block:: python
  
  def accumulate(F, n_divide):
      """Accumulate the total force on the inner fixed disk."""
      F += np.sum(p.get_boundary_forces("Inner"), axis=0) / n_divide
  
  

Time integration loop
---------------------

.. code-block:: python
  
  while t < tEnd:
      print(f"{i:4d}, {t:.6g}/{tEnd:.6g}")
      # Add new particles if the column becomes too low
      nonzero = p.r()[:, 0] != 0
      position = p.position()[nonzero, :]
      if np.amax(position[:, 1]) + r < 1:
          for xi, yi in zip(x, y):
              p.add_particle((xi, yi), r, np.pi * r**2 * rhop)
      # Remove particles below the lower limit
      p.remove_bodies_flag(p.body_position()[:, 1] > -0.6)
      # Write outputs
      if i % outf == 0:
          p.write_mig(outputdir, t)
      # Compute contact and body forces
      mass = np.pi * p.r() ** 2 * rhop
      F = np.zeros(2)
      time_integration.iterate(
          None,
          p,
          dt,
          1,
          contact_tol=1e-3 * r,
          external_particles_forces=g * mass,
          after_sub_iter=lambda n_divide: accumulate(F, n_divide),
      )
      # Write drag results
      with open(csv_file, "a") as file1:
          file1.write(f"{F[0]};{F[1]};{t}\n")
      # Advance time
      t += dt
      i += 1
  

Plot
----
.. code-block:: shell

 python3 -m migflow.plot.migplot output --actors particles

