**Harper is a manufacturer of vertical free-fall reactors for a large number of materials. As temperature increases, CAPEX becomes more significant. As a result it is necessary to use a modeling tool to correctly size the equipment for the desired production rate.**

This article discusses the model, paying attention to the rigorous algorithms and the simplifying assumptions. The model takes into consideration, semi-opaque clouds (efficiency vs. efficacy), heavy particles (acceleration), families (in case sticking is anticipated), collisions between particles in varied size (or shape) and temperature dependent characteristics of the solids and gases (to account for phase changes).

The aim is to form a physically rigorous model in order simulate powder heating with a range of particle sizes. It is possible to use this as a design tool to estimate a system’s productivity.

## Limits

The assumptions made are as follows:

- The model uses a uniform particle distribution across the tube diameter
- A uniform axial gas velocity is also assumed by model, though this is in reality not the case
- Modeling the impact of falling particles on a developing velocity profile was considered highly complex at this stage
- Axial radiation is neglected in the radiation heat transfer calculation. This is not reasonable in the case of large diameter tubes
- The problem of heat transfer stops at the surface of the particle
- The relaxation time for temperature gradients R
^{2}/α may be on an identical order as the free fall time. In those cases the model must be extended along the radius of the particles

## Model

It is possible to model the problem physics to specific limits. The model is in the form of a explicit numerical solution of a system of nonlinear partial differential equations. In the case of explicit integration, any variable’s future value is a function of the present state.

While reducing the time step, the accuracy is enhanced. For specific variables it is easier to integrate directly over space steps in the place of time.

## The Force Balance

Weight = (mass of particle - mass of displaced fluid) x g (gravitational acceleration constant, 9.806 m/s^{2}).

The other force is a drag from friction and pressure between the fluid(gas) and the particle because of the different velocity of the particles and the gas. The drag coefficient (C_{D}) is a dimensionless number enabling the force and the Reynolds number to be related to each other.

Re is also dimensionless, which is the ratio of inertial to viscous forces exerted on the particle by the flow surrounding the particle. C_{D} and the drag force, F_{D} are related and is in the contrary direction of the relative velocity.

F_{D}= C_{D} A_{P} ρ v^{2} / 2

Where A_{P} is the projected area of the particle in the direction of movement

A_{P}=πR^{2} for a sphere, ρ is the density of the fluid or in certain cases the fluid film surrounding the particle and v is the relative velocity between the fluid and the particle.

The Reynolds number is Re = ρ v D / μ. Where D is the diameter of the particle and μ is the fluid’s dynamic viscosity of the fluid, in certain cases the fluid film, around the particle.

Wilke’s mixing rule is used for viscosity

μ_{mixture} = Σ(y_{i}μ_{i}/Σ_{j}Φ_{ij})

Φ_{ij, μ} = 8^{-½} (1+ M_{i}/M_{j})^{-½} (1 + (μ_{i}/μ_{j})^{½} / (M_{i}/M_{j}) -^{¼}))^{2}

Where M_{i} is the molecular weight of species i. By volume, the example mixture is 80%Ar, 20% H_{2}.

The raw material may be crystalline, crushed or spherical. The sphericity Ψ of a particle is the ratio of a sphere’s surface area with the same volume as the given particle to the particle’s surface area.

By definition, the sphericity of a sphere is 1. The shapes such as octahedra, tetrahedral and cubes that one may anticipate to find in a crushed hard material are 0.846, 0.671 and 0.806 respectively.

The sphericity is a parameter in the association between C_{D} and Re.

C_{D} = 24/Re (1 + 8.1716 e^{-4.0655 Ψ} Re^{0.0964 + 0.5565 Ψ} ) + 73.69 e^{- 5.0748 Ψ }Re / (Re + 5.378 e^{6.2122 Ψ})

In the case of spherical particles the above equation is substituted by the much simpler relation:

C_{D} = 24/Re (1 + 0.14Re^{0.7}) valid for 0.1 < Re<1000

The sphericity changes for simulations extending to melting the particles

The acceleration, a_{i}, is (F_{Di} + (6π(ρ _{i}-ρ _{gas})Di^{3})g)/( 6πρ _{i}D_{i}^{3}). F_{Di}, a_{i}, and g are vectors (F up and g down).

Buoyancy modifies the weight force. Velocity is given by :

v_{i} (y+Δy) = v_{i} (y) + a_{i}Δ _{t}. Substitute Δ_{t} = Δ_{y}/v_{i}.

For v_{i} approaching zero, this is unstable. v_{i}=0 is a plugged reactor.

In order to study the impact of crowding other than the effect it may have on radiation heat transfer, collisions are counted.

The assumption is that all particles in a size-group move at the same speed and reach same temperature along the length and continuous particle size distribution is discretized.

**Figure 1. **Viscosity of Argon H2 Mixture as function of Temperature.

*Image Credit: Harper*

The next assumption which may not be true is a uniform gas velocity across the tube cross-section and the velocity of the particles in a particular size range is the same as all other particles in the group at that position.

However, different group particles may move at different velocities.

For example, a large particle group i moving with respect to a large number of smaller group j particles on an average will collide with the volumetric number density of j particles times the volume swept out by the i particle in a specific time increment Δ_{t}. The length increment is Δ_{y}. Hence the time step is Δ_{y}/v_{i} (different for each particle group).

With regards to the collision with j group particles, the volume swept by the i particles are:

V_{c,ij} = π/4(D_{i} + D_{j} ) ^{2}Δ _{t}|v_{j}- V_{i}| = π/4(D_{i} + D_{j})^{2} |v_{j} - v_{i}| Δy/ v_{i}.

**Figure 2. **Continuous Particle Size Distribution.

*Image Credit: Harper*

**Figure 3a. **Centerline Cross Section.

*Image Credit: Harper*

**Figure 3b. **Radial Cross Section.

*Image Credit: Harper*

The volumetric number density of particles of Group j, η _{j} is the number, Δ_{t}N_{j} divided by the volume particles move in a time step Δ_{t}.

η _{j} = Δt m X_{j} / (Δt v_{j} π/4 D_{t}^{2}) = N_{j} / (v_{j} π/4 D_{t}^{2}) (1)

N_{j} = (m X_{i} / (6πρ _{j}D_{j}^{3}))

Where m is the total mass rate (production rate),

X_{j} is the fraction of m in group j, 6πρ _{i}D_{j}^{3} is the particle mass of the j group

Δ _{t} N_{j} is the number of j particles that move across any horizontal plane in the time Δ_{t} the volume of the tube over length Δ_{t} v_{i} is (Δ_{t} v_{i} π D_{t}^{2} / 4) is.

The average collision number that single i particles have with j group particles along length increment Δy is, therefore,

V_{C, ij} x η _{j} = π/4(D_{i}+D_{j})^{2}|v_{j}- v_{i}| Δy/ v_{i} N_{j} / (v_{j} π/4 D_{t}^{2}) = N_{j} (D_{i} +D_{j})^{2}|v_{j}-v_{i}|/(D_{t}^{2} v_{j}v_{i}) Δy

Summing over j, the total number of collisions for one i particle over length Δy, N_{C,i} is

N_{C,i} = {C_{ij}} {N_{j}} Δy,

where C_{ij} = (D_{i} + D_{j})^{2}|v_{j}-v_{i}|/(D_{t}^{2} v_{j}v_{i}) and {N_{j}}= {N_{i}}^{T}.={0, N2, N3,...}, index 0 is the tube.

Thermal energy balance is also explicitly integrated over time.

## Black Bodies

Calculation of radiation heat transfer is done between the furnace and the particles and radiation between particles that because of their different sizes may be at varied temperatures along the free-fall length.

In the case of cylindrical coordinates the areas dA are simply R_{t}d αdz.

In the case of empty tube the integrand is dA_{1} cos(θ_{1}) dA_{2} cos (θ_{2})/pr_{12}^{2}. The cos(θ_{2}) term is the projection of dA_{2} in the direction of dA_{1}. It is possible to prove that cos(θ_{i}) = cos(θ_{2}). The vector from dA_{1} to dA_{2} is r_{12}.

Based on the number density η _{i}, the volumetric density β is defined:

β = Σ _{i} β_{i} = Σ_{i} η_{i} π D_{i}^{2} / 4, β has units of L^{-1} (=L^{2}/L^{3}).

The variation in the intensity I of collimated light is dI = -β I dx.

Integrating over a finite length L, I(L) = I(0) e^{-βL}. A new integrand is obtained by applying the shading of the particles to the radiation view factor calculation

exp(-β|r_{12}|)(|r_{12}•R_{1}|/ |r_{12}| |R_{1}|)² (1/πr_{12}²) R_{t}dα_{1}dz_{1} R_{t}dα _{2}dz_{2}.

A Monte Carlo style algorithm was integrated to this and the result was verified against published values. The points are generated using r/R intercepts from the referenced graph. The calculation of the curve is done with b=0 at varied L/R:

A different cloud density measure is used which is related to β by this simple relation:

A _{tube}/A _{total} = π D_{t} dy / (π D_{t} dy + π/4 D_{t}^{2} dy x 4β) = 1/(1 + βD_{t}).

A little manipulation helps obtain a relationship where f(F_{tt}) and g(A _{tube}/A _{total}) satisfy f(1) =g(1) = 0.

**Figure 4. **Cylindrical Geometry for Integration.

*Image Credit: Harper*

By definition, all views from any giving body sum up to one.

Σ _{j} F_{ij} = 1 (2)

The view which body 1 has of body 2, F_{12} is associated with the view body 2 has of body 1.

A_{1}F_{12} = A_{2}F_{21} or F_{21} = (A_{1}/A_{2}) F_{12} (3)

For determining heat transfer, each group by particle size is treated as one body. The assumption is that particle distribution is uniform across the cross section of the tube. Based on their areas, the view from the tube to the particles is distributed among the particle sets The total particle area is Σ _{j≠i} A_{j}.

The area fraction of a single group, i, is A_{i}/Σ _{j≠i} A_{j}.

**Figure 5. **View Factor of Tube on Tube.

*Image Credit: Harper*

**Figure 6. **View Factor as a Function of Tube Diameter.

*Image Credit: Harper*

**Figure 7. **Linear Fit of View Factor to the Fractional Area.

*Image Credit: Harper*

F_{tp} = 1 – F_{tt}

from equation 2. To simplify the subscripts let F _{ti} = F _{tp}, i,

A_{t} F _{ti} = A _{t}F _{tp} A_{i}/Σ _{j≠1} A_{j},

by distribution over areas argument

F_{pt} = F_{tp} A_{t}/Σ _{j} A_{j},

from equation 3

A_{i} F_{it} = A_{t} F _{Ti} = A _{t}F _{tp} A_{i}/Σ _{j≠1} A_{j},

F_{pp} = 1 – F_{pt}

from equation 2

A_{i} F_{ij} = A_{i}F_{pp} A_{j}/Σ _{j≠1} A_{k} ,

note F_{ii} ≠ 0, group i particles see each other.

Highly dense clouds of particles βD_{t}>>1 may require this work to be extended to verify for radial gradient. The radial gradient, however was not used in the results.

**Figure 8a. **View Factor of Tube to the Tube's Core as a Function of Fractional Area.

*Image Credit: Harper*

## Grey Bodies

The above formulae are for black particles (emissivity = 1). For grey particles (emissivity: independent of wavelength but less than 1) it is important to account for reflected light.

The radiosity is given by:

*Image Credit: Harper*

11 with an over bar is Fractional Area shorthand for A_{1}F_{11}, ε_{1} is the emissivity of body 1 and the reflectivity is ρ1 (ρ = 1-ε, yet another ρ). E_{i} is the blackbody radiative flux from i, σT^{4} where σ is the Stefan-Boltzmann constant.

The system of equations is solved for the W_{i} and with this the net radiation to the surface i is obtained:

The matrix and vectors can be constructed using index vectors, i = {1, 2, ...} and j = i^{T}

{W_{i}}^{T} = {Ω_{ij}}^{-1} {Δ_{i}}^{T}

Ω_{ij} = A_{j}B_{i}F_{ij} - δ_{ij} A_{i}/(1- &epsilon_{i})

Δ_{i} = -A_{i} ε_{i}/(1-ε_{i}) T_{i}^{4}

A_{j} = A_{i}^{T}, A_{i} = {A _{Tube}, A _{Particle group 1}, A _{Particle group 2…}}

B_{i} = {B_{i}}: B_{i=1} = 1, B_{i>1} = A_{i}/ (Σ_{k≠1}A_{k})) …

F_{ij} = {F_{ij}}: F_{i=1, j=1}= F_{TT}, F_{i>1, j=1}= FTP, F_{i=1, j>1}= F_{PT}, F_{i>1, j>1}= F_{PP}

δ_{ij} = 1 if i=j, 0 otherwise

Notes:

The “i” and “j” here include the tube (i, j = 1).

The W in the reference and here differ by a factor of σ.

The radiation heat transfer to a single particle is Q_{rad, i} = σπD_{i}^{2}ε_{i}/(1- ε_{i}) (W_{i} – T_{i}^{4})

**Figure 8b. **View Factor of Shell to the Tube's Core as a Function of Fractional Area.

*Image Credit: Harper*

Flow-dependent empirical correlations are used for calculating the convective heat transfer between the fluid and the particles. It is assumed that the flow is undisturbed by other particles Most simulated cases have a high volumetric loading of below a few percent. It is also assumed that the gas is non-radiating. For systems where gas is flowing, it is important to track or model the gas temperature along the length also.

In case the velocity of gas is zero there cannot be any net heat transfer to or from the gas. This is mathematically stated as the "the sum of the heats transferred to the gas is zero" or: Σ_{i} h_{i}A_{i} (T_{i}-T_{g}) = 0, i h_{i}, A_{i}, and T_{i} are the local heat transfer coefficient, area and temperature of surface "i". Note, again: "i" includes the tube.

The no-net-flow gas temperature is

T_{g, NF} = Σ_{i} (h_{i}A_{i}T_{i})/ Σ_{j} (h_{i}A_{i}).

The gas allows heat transfer between different surfaces in parallel with radiative heat transfer. The process can use flowing gas in different ways. The residence time in the hot zone is extended by gas flowing up past the falling particles. Gas flowing down along with the particles can bring down the velocity range and the number of collisions.

For flowing gas, an explicit integration is more complex. The solution method which works for both extremes and the range between them is to use an integrated form:

m. _{g} C_{p,g} dT_{g} = Σ_{i} h_{i}dA_{i} (T_{i}-T_{g})

Where dA_{i} = (A_{i}/Δy)dz, m. _{g} and C_{p,g} are the mass rate and heat capacity of the gas, respectively. This can be integrated over 0

T_{g}(y+Δy) = T_{g}(y) + (T_{s} - T_{g}(y))(1-exp(-HA/ m._{g} C_{p,g}))

HA = Σ_{i}(h_{i}A_{i}), T_{s} = S_{i} (h_{i}A_{i}T_{i})/HA.

The results do not show sensitivity to if "i" refers to the local temperature or upstream. It is more important to track heat carried by gas for high flow cases.

The correlation used to get the gas-powder heat transfer coefficient is:

Nu_{i} = 2 + 0.6Re_{i}^{½} Pr^{1/3}

Re is previously defined. The Prandtl number Pr ( = C_{p} μ/k) - is considered relatively constant since the thermal conductivity and viscosity of a gas show similar temperature dependence. In the example, the difference over 3000K was around 3%.

h_{i}=Nu_{i}k_{i} /D_{i}, k_{i} = k_{T0} (½/ (T_{i}+T_{g})/T_{o})^{kA}

For cases wherein A_{particies}>A_{tube} and h_{shape} = Nu k_{gas}/D_{shape} yielding h_{particle}>>h_{tube}, it is possible to neglect convection to and from the tube wall.

The particle temperature as a function of time and/or distance is, in explicit numerical integration, a simple algebraic relationship:

**Figure 9. **Velocity and Particle Count for various Particle Sizes.

*Image Credit: Harper*

T(t+Δt) = T(t) + Δt(net heat transfer rate to particle)/(particle mass x heat capacity)

This model allows the passage of materials through energetic phase transitions. This is mathematically captured by making the heat capacity a function of temperature. An enthalpy plot versus temperature is fit with a C_{p} (T). The iteration time steps need to be very small and phase temperature band wide enough so that as the element passes through the phase change, the temperature changes.

## Results

The model helps obtain graphs showing all input parameters. For example Figure 10 shows along with the temperature for each size group and the boundary conditions, the material thermodynamics, the particle size distribution, the particle velocities for each size, the collision count and the full set of input parameters.

**Figure 10**

*Image Credit: Harper*

## About Harper International

Harper International is a global leader in complete thermal processing solutions and technical services essential for the production of advanced materials. From concept to commercialization, from research scale to full production line operations, Harper is perpetually on the cutting edge of the most innovative furnace and oven designs in the world. For decades, we have pioneered some of the most unique, customized systems available, with a focus on processing materials at high temperatures up to 3000°C and in non-ambient atmospheres.

Harper serves advanced, cutting-edge material markets including Fibers & Filaments, Powders, Metal Oxides, Technical Ceramics, Rare Earths, Graphene, Energy Device Materials and Nuclear Materials. Our support to these emerging industries begins in early stages of research and development, whether at corporate R&D centers, universities, government institutions, or start- ups. Harper is a partner through the entire development process assisting in the scale up and commercialization of advanced materials that will change our everyday lives.

This information has been sourced, reviewed and adapted from materials provided by Harper International.

For more information on this source, please visit Harper International.