Download :download:`this testcase <3d_depot_grid.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
  

Tridimensional Particle Sedimentation
=====================================
This example illustrates the sedimentation of dense solid particles in a
viscous fluid under gravity in three dimensions using MigFlow.
Particles are initialized in a grid packing.

Keywords
--------
DEM, FEM, Generation

Description
-----------
The test demonstrates:
- How to generate a 3D rectangular fluid domain with named boundaries in Gmsh.
- How to couple the particle solver (`scontact`) and the fluid solver (`fluid`).
- How to export derived fields (pressure, velocity, porosity, dynamic pressure) for post-processing in Paraview or similar tools.

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

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, 0])  # gravity
  r = 2e-3  # particles radius
  rho = 1000  # fluid density
  rhop = 1500  # particles density
  nu = 1e-6  # kinematic viscosity
  Atwood = (rhop - rho) / (rhop + rho)
  print(f"Atwood number: {Atwood:1.4f}")
  

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

.. code-block:: python
  
  w = 0.4  # domain width
  h = 0.6  # domain height
  d = 10 * r  # domain depth
  lx = w  # Particles width
  ly = h / 6  # Particles height
  lz = d  # Particles depth
  lc = d / 2  # mesh size
  

Mesh generation
---------------

.. code-block:: python
  
  gmsh.initialize()
  gmsh.model.add("box3d")
  gmsh.model.occ.add_box(-w / 2, 0, -d / 2, w, h, d)
  gmsh.model.occ.synchronize()
  
  
  def get_surface(x0, x1, y0, y1, z0, z1, eps=1e-6):
      r = gmsh.model.get_entities_in_bounding_box(
          x0 - eps, y0 - eps, z0 - eps, x1 + eps, y1 + eps, z1 + eps, 2
      )
      return [tag for dim, tag in r]
  
  
  gmsh.model.add_physical_group(
      2, get_surface(-w / 2, +w / 2, h, h, -d / 2, +d / 2), name="Top"
  )
  gmsh.model.add_physical_group(
      2, get_surface(-w / 2, +w / 2, 0, 0, -d / 2, +d / 2), name="Bottom"
  )
  gmsh.model.add_physical_group(
      2, get_surface(-w / 2, -w / 2, 0, h, -d / 2, +d / 2), name="Left"
  )
  gmsh.model.add_physical_group(
      2, get_surface(+w / 2, +w / 2, 0, h, -d / 2, +d / 2), name="Right"
  )
  gmsh.model.add_physical_group(
      2, get_surface(-w / 2, +w / 2, 0, h, +d / 2, +d / 2), name="Front"
  )
  gmsh.model.add_physical_group(
      2, get_surface(-w / 2, +w / 2, 0, h, -d / 2, -d / 2), name="Rear"
  )
  gmsh.model.add_physical_group(3, [1], name="domain")
  
  # --- Uniform mesh control ---
  gmsh.model.mesh.set_algorithm(3, 1, 10)
  gmsh.option.setNumber("Mesh.CharacteristicLengthMin", lc)
  gmsh.option.setNumber("Mesh.CharacteristicLengthMax", lc)
  gmsh.option.setNumber("Mesh.CharacteristicLengthFromCurvature", 0)
  gmsh.option.setNumber("Mesh.CharacteristicLengthFromPoints", 0)
  gmsh.option.setNumber("Mesh.CharacteristicLengthExtendFromBoundary", 0)
  gmsh.model.mesh.generate(3)
  
  

Function to generate a rectangular grid of particles
----------------------------------------------------

.. code-block:: python
  
  def gen_rect3d(origin, w, h, d, step):
      """Generate coordinates on a rectangular grid
      Keyword arguments:
          origin -- origin of top left corner
          w : width
          h : height
          d : depth
          step : maximal radius
      """
      eps = 1e-8
      x = np.arange(-w / 2 + step - eps, w / 2 - step + eps, 2 * step) + origin[0]
      y = np.arange(origin[1] - step - eps, origin[1] - h + step + eps, -2 * step)
      z = np.arange(-d / 2 + step - eps, d / 2 - step + eps, 2 * step) + origin[2]
      x, y, z = np.meshgrid(x, y, z)
      return x.reshape(-1), y.reshape(-1), z.reshape(-1)
  
  

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

.. code-block:: python
  
  p = scontact.ParticleProblem(3)
  p.load_msh_boundaries(None)
  x, y, z = gen_rect3d([0, h, 0], lx, ly, lz, r)
  for xi, yi, zi in zip(x, y, z):
      p.add_particle((xi, yi, zi), r, 4 / 3 * r**3 * np.pi * rhop)
  

Fluid Problem Initialization

.. code-block:: python
  
  f = fluid.FluidProblem(3, g, nu * rho, rho)
  f.load_msh(None)
  f.set_wall_boundary(["Left", "Right", "Front", "Rear", "Top", "Bottom"])
  f.set_mean_pressure(0)
  
  

Function to get derived fields for output
-----------------------------------------

.. code-block:: python
  
  def get_fields(fluid):
      """Return derived output fields for visualization."""
      y = fluid.coordinates_fields()[fluid.pressure_index()][:, 1].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[1] * (y - h / 2), p1_element),
      }
  
  

Simulation Parameters
---------------------

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

Simulation Loop
---------------
Time integration of coupled fluid–particle motion.

.. code-block:: python
  
  mass = 4 / 3 * np.pi * p.r() ** 3 * rhop
  while t < tEnd:
      print(f"t/tEnd={t:1.4f}/{tEnd:1.4f}, i={i:1d}")
      f.set_particles(
          p.delassus(),
          p.volume(),
          p.position(),
          p.velocity(),
          p.omega(),
          p.contact_forces(),
      )
      if i % outf == 0:
          p.write_mig(outputdir, t)
          f.write_mig(outputdir, t, get_fields(f))
      f.implicit_euler(dt)
      fext = g * mass + f.compute_node_force()
      time_integration._advance_particles(p, fext, dt, 1, 1e-3 * r)
      i += 1
      t += dt
