Free-Fall Heating of Particle to High Temperatures

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²/α 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²).

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

Where A_P is the projected area of the particle in the direction of movement

A_P=πR² 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) -^¼))²

Where M_i is the molecular weight of species i. By volume, the example mixture is 80%Ar, 20% H₂.

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³)g)/( 6πρ _iD_i³). 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 ) ²Δ _t|v_j- V_i| = π/4(D_i + D_j)² |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, Δ_tN_j divided by the volume particles move in a time step Δ_t.

   η _j = Δt m X_j / (Δt v_j π/4 D_t²) = N_j / (v_j π/4 D_t²)         (1)

   N_j = (m X_i / (6πρ _jD_j³))

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

X_j is the fraction of m in group j, 6πρ _iD_j³ 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² / 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)²|v_j- v_i| Δy/ v_i N_j / (v_j π/4 D_t²) = N_j (D_i +D_j)²|v_j-v_i|/(D_t² v_jv_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)²|v_j-v_i|/(D_t² v_jv_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_td αdz.

In the case of empty tube the integrand is dA₁ cos(θ₁) dA₂ cos (θ₂)/pr₁₂². The cos(θ₂) term is the projection of dA₂ in the direction of dA₁. It is possible to prove that cos(θ_i) = cos(θ₂). The vector from dA₁ to dA₂ is r₁₂.

Based on the number density η _i, the volumetric density β is defined:

   β = Σ _i β_i = Σ_i η_i π D_i² / 4, β has units of L^-1 (=L²/L³).

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₁₂|)(|r₁₂•R₁|/ |r₁₂| |R₁|)² (1/πr₁₂²) R_tdα₁dz₁ R_tdα ₂dz₂.

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² 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₁₂ is associated with the view body 2 has of body 1.

   A₁F₁₂ = A₂F₂₁ or F₂₁ = (A₁/A₂) F₁₂         (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 _tF _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 _tF _tp A_i/Σ _j≠1 A_j,

   F_pp = 1 – F_pt

from equation 2

   A_i F_ij = A_iF_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₁F₁₁, ε₁ is the emissivity of body 1 and the reflectivity is ρ1 (ρ = 1-ε, yet another ρ). E_i is the blackbody radiative flux from i, σT⁴ 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_jB_iF_ij - δ_ij A_i/(1- &epsilon_i)

   Δ_i = -A_i ε_i/(1-ε_i) T_i⁴

   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≠1A_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²ε_i/(1- ε_i) (W_i – T_i⁴)

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_iA_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_iA_iT_i)/ Σ_j (h_iA_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_idA_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_iA_i), T_s = S_i (h_iA_iT_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_ik_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.