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

Drag on a Fixed Disk in a 2D Immersed Granular Flow
===================================================
This example studies the drag force acting on a fixed circular obstacle
(a “disk”) immersed in a granular flow saturated by a viscous fluid.


Keywords
--------
Drag, DEM, FEM, Boundary Force, 2D-3D Ratio


Description
-----------
The test demonstrates:
- How to set up a mixed granular–fluid simulation in a 2D geometry.
- How to impose an inflow boundary condition to simulate a stream.
- How to insert, remove, and constrain particles dynamically.
- How to compute and record the total drag force acting on a solid obstacle.
A 2D–3D scaling factor is applied to account for the reduced dimensionality.

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

Output Directory and Mesh Generation
------------------------------------
The output directory is created, and a 2D annular mesh is generated using Gmsh.

.. 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
  
  h_p = 0.3  # particle bed height [m]
  w_p = 0.3  # particle bed width [m]
  w = 0.3  # particle bed width [m]
  h = 1.4  # total domain height [m]
  r_inner = 0.025  # radius of the fixed disk [m]
  ratio_2d_3d = 0.7  # 2D-3D correction
  

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

.. code-block:: python
  
  v_imposed = -0.05  # stream velocity
  rhop = 2500  # particle density
  rho = 1000  # fluid density
  r = 3e-3  # particle radius
  g = np.array([0, -9.81])  # gravity
  nu = 1e-6  # fluid kinematic viscosity
  mu = rho * nu
  friction = 0.2
  
  

Particle Problem Initialization
-------------------------------
Load mesh boundaries and define contact properties.

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

Particle Generation
-------------------
The following function generates particle coordinates on a regular grid
over a rectangular region.

.. code-block:: python
  
  def gen_rect(origin, w, h, step):
      """Generate coordinates on a rectangular grid
      Keyword arguments:
          origin -- origin of top left corner
          w : width
          h : height
          step : maximal radius
      """
      eps = 1e-4 * step
      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.reshape(-1), y.reshape(-1)
  
  
  x, y = gen_rect([0, -h / 2], w, 4 * h_p, r)
  for xi, yi in zip(x, y):
      if xi**2 + yi**2 > (r_inner + r) ** 2:
          p.add_particle((xi, yi), r, r**2 * np.pi * rhop, "Sand")
  x, y = gen_rect([0, 0.5 * h / 1.4], w, 0.5 * h_p, 1.2 * r)
  for xi, yi in zip(x, y):
      p.add_particle((xi, yi), r, r**2 * np.pi * rhop, "Sand")
  
  

Fluid Problem Initialization
----------------------------
The fluid problem is loaded from the same mesh and coupled to the particle phase.

.. code-block:: python
  
  f = fluid.FluidProblem(2, g, nu * rho, rho)
  f.load_msh(mesh_filename)
  f.set_wall_boundary("Inner")
  f.set_wall_boundary("Left", velocity=[0, v_imposed])
  f.set_wall_boundary("Right", velocity=[0, v_imposed])
  f.set_open_boundary("Bottom", velocity=[0, v_imposed])
  f.set_strong_boundary("Left", velocity=[0, v_imposed])
  f.set_strong_boundary("Right", velocity=[0, v_imposed])
  f.set_strong_boundary("Bottom", velocity=[0, v_imposed])
  f.set_open_boundary("Top", pressure=0)
  

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

.. code-block:: python
  
  dt = 5e-3  # time step
  t = 0  # initial time
  tEnd = 1.5  # final time
  i = 0  # iteration number
  outf = 2  # iterations between data frames
  
  

Computation Loop
----------------
Function to accumulate the force applied to the inner disk by the contact networks.

.. code-block:: python
  
  def accumulate(F, n_divide):
      F += np.sum(p.get_boundary_forces("Inner"), axis=0) / n_divide
  
  
  while t < tEnd:
      print(f"{i:4d}, {t:.6g}/{tEnd:.6g}")
      # Remove particles that have fallen outside the domain
      alert = 4 * r
      keep = (p.body_position()[:, 1] > -h / 1.4 * 0.55) & (
          p.body_position()[:, 1] < h / 2 - alert
      )
      p.remove_bodies_flag(keep)
  
      # Constrain particles at the bottom
      nonzero = p.r()[:, 0] != 0
      position = p.position()[nonzero, :]
      p.body_invert_mass()[p.body_position()[:, 1] < -0.5] = 0
      p.body_invert_inertia()[p.body_position()[:, 1] < -0.5] = 0
      a = p.body_velocity()[p.n_bodies() - position.shape[0] :, :]
      b = p.body_omega()[p.n_bodies() - position.shape[0] :, :]
      c = p.body_invert_mass()[p.n_bodies() - position.shape[0] :, :]
      a[(c == 0)[:, 0], :] = [0, v_imposed]
      b[(c == 0)[:, 0], 0] = 0
  
      # Insert new particles at the top when needed
      if np.amax(position[:, 1]) + r < 0.5:
          for xi, yi in zip(x, y):
              p.add_particle((xi, yi), r, r**2 * np.pi * rhop, "Sand")
  
      # Set particles to the FEM module
      f.set_particles(
          p.delassus(),
          p.volume() * ratio_2d_3d,
          p.position(),
          p.velocity(),
          p.omega(),
          p.contact_forces(),
      )
  
      # Write output files
      if i % outf == 0:
          p.write_mig(outputdir, t)
          f.write_mig(outputdir, t)
  
      # Iterate
      F = np.zeros(2)  # to store inner disk forces
      mass = np.pi * p.r() ** 2 * rhop
      time_integration.iterate(
          f,
          p,
          dt,
          10,
          contact_tol=1e-3 * r,
          external_particles_forces=g * mass,
          after_sub_iter=lambda n_divide: accumulate(F, n_divide),
          use_predictor_corrector=False,
      )
  
      with open(csv_file, "a") as file1:
          F += f.boundary_force()[0, :]
          file1.write(str(F[0]) + ";" + str(F[1]) + ";" + str(t) + "\n")
      t += dt
      i += 1
  

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

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

