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

Darcy Flow Through a Bed of Circular Discs
==========================================
This example simulates Darcy flow through a packed bed of circular
disc-shaped particles in two dimensions. A pressure gradient is applied
across the domain and the resulting flow field through the porous medium
is computed.

Keywords
--------
FEM, DEM, Darcy, Porous Media

.. code-block:: python
  
  import os
  import gmsh
  import shutil
  import numpy as np
  from migflow import fluid, scontact
  
  outputdir = "output"
  shutil.rmtree(outputdir, ignore_errors=True)
  os.makedirs(outputdir)
  

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

.. code-block:: python
  
  r = 5e-3  # grains radius
  L = 25 * 2 * r  # domain width
  H = 25 * 2 * r  # domain height
  mesh_size = r / 4  # 0.1                             # mesh size
  # mesh_size = 0.1                             # mesh size
  # mesh_size = L/4                             # mesh size
  origin = np.array([-L / 2, -L / 2])  # mesh origin
  

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

.. code-block:: python
  
  rho = 1000  # fluid density
  rhop = 1500  # grains density
  mu = 1  # dynamic viscosity
  compacity = 0.5  # Compactness
  porosity = 1 - compacity  # porosity
  p1 = 0
  p0 = 10.0
  g = np.array([-(p1 - p0) / (rho * L), 0])  # gravity
  
  
  gmsh.initialize()
  gmsh.model.add("mesh")
  
  x0 = origin
  x1 = origin + np.array([L, 0])
  x2 = origin + np.array([L, H])
  x3 = origin + np.array([0, H])
  x = np.array([x0, x1, x2, x3])
  edges = np.array([[0, 1], [1, 2], [2, 3], [3, 0]])
  for xi in x:
      gmsh.model.geo.add_point(xi[0], xi[1], 0, mesh_size)
  for edge in edges:
      gmsh.model.geo.add_line(edge[0] + 1, edge[1] + 1)
  gmsh.model.geo.add_curve_loop([1, 2, 3, 4], 1)
  gmsh.model.geo.add_plane_surface([1], 1)
  gmsh.model.geo.synchronize()
  
  
  def get_line(x0, x1, eps=1e-6):
      bbentities = 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 bbentities)
  
  
  # gmsh.model.add_physical_group(1,get_line([origin[0],origin[1]], [origin[0]+L,origin[1]]), name="Bottom")
  # gmsh.model.add_physical_group(1,get_line([origin[0],origin[1]+L], [origin[0]+L,origin[1]+L]), name="Top")
  gmsh.model.add_physical_group(
      1, get_line([origin[0], origin[1]], [origin[0], origin[1] + L]), name="Left"
  )
  gmsh.model.add_physical_group(
      1,
      get_line([origin[0] + L, origin[1]], [origin[0] + L, origin[1] + L]),
      name="Right",
  )
  gmsh.model.add_physical_group(2, [1], name="domain")
  transform_y = np.array([[1, 0, 0, 0], [0, 1, 0, -H], [0, 0, 1, 0], [0, 0, 0, 1]])
  transform_x = np.array([[1, 0, 0, +L], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]])
  
  gmsh.model.mesh.set_periodic(1, [1], [3], transform_y.flatten().tolist())
  gmsh.model.mesh.set_periodic(1, [2], [4], transform_x.flatten().tolist())
  gmsh.model.mesh.set_size_callback(lambda dim, tag, x, y, z, lc: mesh_size)
  gmsh.model.mesh.generate(2)
  p = scontact.ParticleProblem(2)
  p.load_msh_boundaries(None)
  
  # Definition of the points where the grains are located
  # The particles are first placed on a regular grid
  # Then they are centered in the domain
  a_g = np.pi * r**2
  a_b = L**2
  N = max(compacity * a_b / a_g, 1)
  
  e = L / N**0.5
  x = np.linspace(-L / 2, L / 2, int(N**0.5))
  y = np.linspace(-H / 2, H / 2, int(N**0.5))
  # x = np.arange(-L/2, L/2-r, e)
  # y = np.arange(-L/2, L/2-r, e)
  x, y = np.meshgrid(x, y)
  for xi, yi in zip(x.flat, y.flat):
      p.add_particle((xi, yi), r, r**2 * np.pi * rhop)
  
  f = fluid.FluidProblem(2, g, mu, rho)
  f.load_msh(None)
  f.set_mean_pressure(0)
  # f.interpolate(velocity_x=U)
  
  

Numerical parameters
--------------------

.. code-block:: python
  
  dt = 1e-4  # time step
  tEnd = 100 * dt  # final time
  outf = 10  # number of iterations between output files
  t = 0.0
  i = 0
  

Time Integration
----------------

.. code-block:: python
  
  
  
  def get_fields(fluid):
      """Return derived output fields for visualization."""
      x = fluid.coordinates_fields()[fluid.pressure_index()][:, 0].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),
          "u_solid": (fluid.u_solid(), p1_element),
          "dynamic_pressure": (fluid.pressure() - rho * g[0] * x, p1_element),
      }
  
  
  while t < tEnd:
      print(f"t/tEnd={t:1.4f}/{tEnd:1.4f}, i={i:1d}")
      if i % outf == 0:
          p.write_mig(outputdir, t)
          f.write_mig(outputdir, t, get_fields(f))
      contacts = f.compute_node_force() if i > 0 else np.zeros_like(p.velocity())
      f.set_particles(
          p.delassus(), p.volume(), p.position(), p.velocity(), p.omega(), -contacts
      )
      f.implicit_euler(dt)
      i += 1
      t += dt
  
  node_volume = f.node_volume()
  
  U = np.sum(f.velocity()[:, 0] * node_volume) / np.sum(node_volume)
  V = U / porosity
  Re = V * rho * 2 * r / mu
  print("Re : %2.g ---- t : %2.g ---- tEnd : %2.g" % (Re, dt, tEnd))
  
  voidage = porosity ** (-1.8)
  cross_section = 2 * r
  Re = U * rho * cross_section / mu
  Cd = voidage * (0.63 + 4.8 * (porosity * Re) ** (-0.5)) ** 2
  gamma = Cd * 0.5 * rho * cross_section * U
  vol = np.pi * r**2
  N = compacity / vol
  Dp_drag = -gamma * N * V
  
  dp = f.fields_gradient()[:, -1]
  dp = np.sum(dp * node_volume) / np.sum(node_volume) - rho * g[0]
  
  error = np.sum(abs(dp - Dp_drag) / Dp_drag)
  print(dp, Dp_drag, error)
  
  print("Pressure gradient from simulation : ", dp)
  print("Pressure gradient from drag model : ", Dp_drag)
  print("Relative error : ", error)
  

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

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

