Download this testcase.

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.

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.

outputdir = "output" if len(sys.argv) < 2 else sys.argv[1]
shutil.rmtree(outputdir, ignore_errors=True)
os.makedirs(outputdir)

Geometrical Parameters

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

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

Particle Problem Initialization

p = scontact.ParticleProblem(2)

Generate the Silo Geometry

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

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

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

i = 0
t = 0
dt = 2.5e-3
tEnd = 4
outf = 10
nsub = 5

Simulation Loop

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

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