*2.2. DEM Method and Heat Transfer Model*

In the present study, the flow and heat transfer of particles are simulated with discrete element method (DEM). For DEM simulations, the Newton's second law [23] is used to calculate the forces on each particle and track the particle motions. For gravity-driven dense particle flow, since the particle filling rate is high and particles move very slowly, the effect of the gas flow between particles is insignificant [24]. Therefore, the gas flow is not considered for the simulations, and the Hertz-Mindlin soft sphere model [25,26] is adopted to simulate the dense particle flow process, where the normal force and tangential force between particles are treated with equivalent spring, damper and slider to simplify the contact between particles. The contact force is calculated according to the normal overlap and tangential displacement of particles. The equations used to calculate the contact force of particle-particle or particle-wall are as follows:

$$F\_{\rm n} = \frac{4}{3} E\_{\rm eq} \sqrt{r\_{\rm eq}} \delta\_{\rm n}^{3/2} + 2 \sqrt{\frac{5}{6}} \frac{\ln \varepsilon}{\sqrt{\ln^2 \varepsilon + \pi^2}} \sqrt{2 E\_{\rm eq} m\_{\rm eq}} \sqrt{r\_{\rm eq} \delta\_{\rm n}} V\_{\rm n}^{\rm rel} \tag{1}$$

$$F\_{\rm t} = \begin{cases} 8G\_{\rm eqq} \delta \sqrt{r\_{\rm eq} \delta\_{\rm n}} + 2\sqrt{\frac{5}{6}} \frac{\ln c}{\sqrt{\ln^2 c + \pi^2}} \sqrt{8G\_{\rm eqq} m\_{\rm eq} \sqrt{r\_{\rm eq} \delta\_{\rm n}}} V\_{\rm t}^{\rm ref} & |F\_{\rm t}| < \mu\_{\rm s} |F\_{\rm n}| \\\ \mu\_{\rm s} F\_{\rm n} \text{sign}(\delta\_{\rm t}) & |F\_{\rm t}| \ge \mu\_{\rm s} |F\_{\rm n}| \end{cases} \tag{2}$$

where *F*n and *F*t are the normal and tangential components of the contact force. *E*eq, *R*eq, and *m*eq are the equivalent Young's modulus, radius, and mass. δ<sup>n</sup> and δ<sup>t</sup> are normal and tangential displacements. *V*n rel and *V*<sup>t</sup> rel are the relative normal and tangential translational velocities. *e* and μ<sup>s</sup> are the restitution coefficient and the translational friction coefficient.

As for heat transfer of particle flow, the heat transfer form particle to particle (p-p) or from particle to the tube wall (p-w) includes heat conduction, thermal convection and thermal radiation, where the heat conduction consists of heat conduction inside particles, physical contact heat conduction and heat conduction through gas film surround particles. In the present study, since particles move quite slowly, the gas convection heat transfer is very small [27]. The thermal convection of gas flow between particles is not considered in the simulations. The heat transfer model adopted in the present paper is based on the following assumptions, (1) the diameter of all particles is the same; (2) the gas heat capacity is ignored [28]; (3) particles are wrapped with gas film and its thickness is 0.1 *d*p [29]; (4) the heat transfer between particles is along the radial direction of particles; (5) the thermal properties of gas and particles are kept constant, as listed in Table 2. The correlation of the heat transfer (*Q*) with temperature difference (Δ*T*), total thermal resistance (*R*total) and time (*t*) is presented in Equation (3), and the total thermal resistance (*R*total) is calculated based on the thermal resistance network, as shown in Figure 2.

$$Q = \Delta t \Delta T / \mathcal{R}\_{\text{total}} \tag{3}$$

**Table 2.** Typical physical properties of gas and particles (Liu et al. [1]).


**Figure 2.** Schematic of thermal resistance network: (**a**) the case for particle–particle; (**b**) the case for particle–wall (*R*1: thermal conduction resistance inside particles; *R*2: physical contact thermal conduction resistance; *R*3: thermal conduction resistance through gas film; *R*4: thermal radiation resistance; "G-C" means contact only with gas film; "P–C" means physical contact).

The thermal conduction resistance includes thermal conduction resistance inside particles (*R*1), physical contact thermal conduction resistance (*R*2) and thermal conduction resistance through gas film (*R*3). According to Fourier's law [30], the thermal resistance *R*<sup>1</sup> and the thermal resistance *R*<sup>3</sup> are formulated as follows:

$$R\_1 = \frac{\sqrt[3]{2} - 1}{2\pi k\_3 r (1 - \cos \alpha)}\tag{4}$$

$$R\_3 = \left[ \int\_{\beta}^{\alpha} \frac{2\pi k\_{\rm g} r^2 \sin\theta \cos\theta}{l\_{\rm ii} - r \sin\theta - \sqrt{r^2 - \left(r \sin\theta\right)^2}} d\theta \right]^{-1} \tag{5}$$

where α is an angle related to the intersection of gas film and β is a starting point angle for the calculation of thermal conduction resistance through gas film, which are both presented in Figure 3 and Equations (6) and (7); *l*ij is the distance for the case of particle-particle or particle-wall, as shown in Figure 3; *r* is the particle radius and *r*<sup>g</sup> is the particle radius plus gas film (*r*<sup>g</sup> = 0.6 *d*p).

$$\alpha = \begin{cases} \arcsin\left(\sqrt{r\_{\text{g}}^2 - \left(\frac{l\_{\text{i}}}{2}\right)^2}/r\right) & (\text{p-p})\\ \arcsin\left(\sqrt{r\_{\text{g}}^2 - l\_{\text{i}}^2}/r\right) & (\text{p-w}) \end{cases} \tag{6}$$

$$\beta = \arccos\left(\frac{l\_{\bar{\imath}\bar{\jmath}} - 6.8 \times 10^{-8}}{2r}\right) \tag{7}$$

when the contact is only with gas film and the gas film overlap is larger than gas free path, the value of β is equal to 0.

The physical contact thermal conduction resistance *R*<sup>2</sup> [31] is formulated as follows:

$$R\_2 = \left[2k\_s \left(\frac{3F\_{\rm n}r}{4E\_{\rm eq}}\right)^{1/3}\right]^{-1} \tag{8}$$

where *F*n is normal direction force and *E*eq is Young's modulus.

According to the Stefan–Boltzmann law, the thermal radiation resistance *(R*4) is presented in Equation (9). *X*ij is the view factor, as presented in Equation (10). σ is the Stefan–Boltzmann constant and ε is the surface emissivity. When the particle contacts with the tube wall, *X*ij is equal to 0.315 [32].

$$R\_4 = \left[\frac{\sigma \left(T\_{\rm i}^2 + T\_{\rm j}^2\right) \left(T\_{\rm i} + T\_{\rm j}\right)}{\left(\frac{1-\varepsilon}{\varepsilon A}\right)\_{\rm i} + \frac{1}{A\_i X\_{\rm ij}} + \left(\frac{1-\varepsilon}{\varepsilon A}\right)\_{\rm j}}\right]^{-1} \tag{9}$$

*Energies* **2020**, *13*, 1961

$$X\_{\overline{\mathbf{q}}} = \begin{cases} \left[1 - \sqrt{1 - \left(\frac{r\_i}{l\_{\overline{\mathbf{q}}}}\right)^2}\right] \left[1 - \sqrt{1 - \left(\frac{r\_i}{l\_{\overline{\mathbf{q}}}}\right)^2}\right] \left(\frac{l\_{\overline{\mathbf{q}}}}{r\_i}\right)^2 & (\text{p-p})\\ 0.315 & (\text{p-w}) \end{cases} \tag{10}$$

At the beginning of DEM simulations, randomly packed high-temperature particles were generated in the particle flow channel. Then, particles start to flow in the channel under gravity. High temperature particles are cooled when they flow around the tube. Meanwhile, new particles were generated and continuously enter into the channel at the inlet. During the whole simulation process, the total particle number inside the channel is kept constant. Typical mechanical parameters for the DEM simulations are listed in Table 3, and the heat transfer coefficient (*h*) of particle flow around the tube is defined as follows:

$$\hbar = \frac{q}{A\_{\text{tube}}(T\_{\text{in}} - T\_{\text{tube}})} \tag{11}$$

where *q* is the heat flux on the tube wall; *A*tube is the heat transfer area of tube wall; *T*in is the particle inlet temperature; *T*tube is the tube wall temperature.

**Figure 3.** Schematic diagram of thermal conduction resistance calculation parameters: (**a**) physical contact for the case of particle-particle and particle-wall; (**b**) contact only with gas film for the case of particle-particle and particle-wall.

**Table 3.** Typical mechanical parameters for discrete element method (DEM) simulations.


The variations of heat flux (*q*) on the tube wall with time for the particle flow around a circular tube at *v*out = 0.5 mm/s are presented in Figure 4. It shows that, when d*q*/d*t* is quite small (d*q*/d*t* < 0.05%), the particle flow and heat transfer around the tube should be quasi-steady. In the present study, the simulation results within the 30 s after the quasi-steady state are extracted for the analysis.

**Figure 4.** Variations of heat flux with time for the particle flow around a circular tube (*v*out = 0.5 mm/s, *d*p = 1.72 mm).
