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

Heat Transfer in a 2D Silo with Granular Flow
=============================================
This example studies heat conduction through a granular material flowing
out of a silo, with thermal interaction between particles.

Keywords
--------
Heat Transfer, DEM, Thermal Conduction, Granular Flow, Silo Discharge

Description
-----------
The test demonstrates:
- How to set up a granular flow simulation in a silo geometry.
- How to model heat transfer and thermal conduction between particles.
- How to initialize particles with temperature fields.
- How to remove particles that exit the system dynamically.
- How to track and record temperature evolution during flow.

.. code-block:: python
  
  
  import os, sys, shutil, time
  import numpy as np
  from migflow import 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)
  

Geometrical Parameters
----------------------

.. code-block:: python
  
  h = 10  # silo height
  w = 20  # silo width
  h_outlet = 4  # silo outlet height
  l_outlet = 2  # silo outlet width
  r = 0.1  # particle radius
  tol = 1e-3 * r
  

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

.. code-block:: python
  
  rho = 1
  cp = 1
  g = np.array([0.0, -9.81])
  

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

.. code-block:: python
  
  p = scontact.ParticleProblem(2)
  

Generate the Silo Geometry

.. code-block:: python
  
  p_tl = np.array([-w / 2, +h / 2])
  p_bl = np.array([-w / 2, -h / 2])
  p_ol = np.array([-l_outlet / 2, -h / 2 - h_outlet])
  p_or = np.array([+l_outlet / 2, -h / 2 - h_outlet])
  p_br = np.array([+w / 2, -h / 2])
  p_tr = np.array([+w / 2, +h / 2])
  
  
  pts = np.array([p_tl, p_bl, p_ol, p_or, p_br, p_tr])
  e0 = np.arange(0, pts.shape[0])
  e1 = np.arange(1, pts.shape[0] + 1) % pts.shape[0]
  masks = np.ones(e0.shape[0], dtype=bool)
  masks[2] = False  # skip outlet edge
  body_bnd = p.add_body([0, 0], 0.0, 0.0)
  for e0, e1, mask in zip(e0, e1, masks):
      if not mask:
          continue
      p.add_segment_to_body(pts[e0], pts[e1], body_bnd, "wall")
  for pt in pts:
      p.add_particle_to_body(pt, 0.0, body_bnd)
  

Particle Initialization via Fast Deposit
----------------------------------------

.. code-block:: python
  
  np.random.seed(0)
  radii = np.random.uniform(0.7, 1, size=5000) * r
  xp, rp = scontact.fastdeposit_2d(pts, radii, "deposit.svg")
  for xi, ri in zip(xp, rp):
      p.add_particle(xi, ri, np.pi * ri**2 * rho)
  
  

Material Properties and Initial Conditions
------------------------------------------

.. code-block:: python
  
  def set_heat_properties(p):
      ks = np.full(p.n_bodies(), 1.0)
      young_modulus = np.full(p.n_bodies(), 1.0)
      poisson_coeff = np.full(p.n_bodies(), 0.3)
      inv_vol = np.zeros(p.n_bodies())
      vol = np.pi * p.r() ** 2
      inv_vol[p.particle_body()] = 1 / vol.reshape(-1)
      inv_vol[body_bnd] = 0.0  # no heat transfer for boundary particles
      return ks, young_modulus, poisson_coeff, inv_vol
  
  
  T0 = 0.0
  tp = np.zeros(p.n_bodies()) + T0
  tp[body_bnd] = 1 + T0
  

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

.. code-block:: python
  
  i = 0
  t = 0
  dt = 2.5e-3
  tEnd = 4
  outf = 10
  nsub = 5

Simulation Loop
---------------

.. code-block:: python
  
  while t < tEnd:
      print(f"{i:4d}, {t:.6g}/{tEnd:.6g}")
  
      keep = np.ones(p.n_bodies(), dtype=bool)
      keep[p.body_position()[:, 1] < -h / 2 - h_outlet - 10 * r] = False
      p.remove_bodies_flag(keep)
      # Remove corresponding temperature entries
      tp = tp[keep]
      ks, young_modulus, poisson_coeff, inv_vol = set_heat_properties(p)
  
      if i % outf == 0:
          p.write_mig(
              outputdir, t, extra_fields={"temperature": tp[p.particle_body()] - T0}
          )
      # Dynamics
      tic = time.process_time()
      mass = rho * p.volume()
      time_integration._advance_particles(p, mass * g, dt, 1, tol)
      print(f"\tNSCD : time ----- {time.process_time()-tic}")
      tic = time.process_time()
      # Heat Transfer
      told = tp.copy()
      for isub in range(nsub):
          tp += (
              dt
              / nsub
              * p.contact_heat(ks, young_modulus, poisson_coeff, tp)
              * inv_vol
              / (rho * cp)
          )
      print(f"\tHEAT : time ----- {time.process_time()-tic} ----- {np.max(tp-told)/dt=}")
      # Update time
      t += dt
      i += 1
  

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

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

