Lattice Boltzmann Simulation of a Gas-to-Solid Reaction and Precipitation Process in a Circular Tube

Matthew D. Lindemer , Suresh G. Advani and Ajay K. Prasad

Abstract: The lattice Boltzmann method (LBM) is used to simulate the growth of a solid-deposit on the walls of a circular tube resulting from a gas-to-solid reaction and precipitation process. This process is of particular interest for the design of reactors for the production of hydrogen by the heterogeneous hydrolysis of steam with Zn vapor in the Zn/ZnO thermochemical cycle. The solid deposit of ZnO product on the tube wall evolves in time according to the temporally- and axially-varying convective-diffusive transport and reaction of Zn vapor with steam on the solid surface. The LBM is wellsuited to solving problems with coupled flow, heat and mass transfer in a time-evolving domain. Here, a D2Q9 axisymmetric multiple-relaxation-time (MRT) lattice Boltzmann scheme is used to simulate incompressible fluid transport while a D2Q5 axisymmetric MRT lattice Boltzmann scheme is used to simulate the convective-diffusive transport of Zn vapor. The model is first validated against several analytical solutions, followed by a parametric study to understand the effect of Reynolds, Schmidt, and Damk?hler numbers on the time evolution of the ZnO deposition profile along the tube axis. The axial location of the fastest deposition is found to increase with increasing Peclet number, and decrease with increasing Damk?hler number, with no independent effect from the Schmidt number.When the reaction kinetics are assumed to increase along the tube axis due to nonisothermal tube wall temperature, a second peak in the deposition profile can be observed for sufficiently low values of Da/Pe.

Keywords: Lattice Boltzmann methods, reactive flow, heterogeneous reaction,precipitation, solar hydrogen production.

Nomenclature

D mass diffusivity

e discrete velocity vector

f fluid distribution function

g concentration distribution function

G flow source term

K kinetic constant

L length of modeled section of tube

M collision matrix

N number of lattice units

P pressure

Q volumetric flow rate

r radial coordinate

rffreaction rate

R tube radius

R mass transfer source term

S source term

t time

u hydrodynamic velocity

w weight function

W wall location function

z axial coordinate

z0beginning of non-isothermal tube section

α discrete velocity index

Γ oscillation time period

δ Kronecker delta

Δ boundary intersection fraction

θ local angle of solid boundary

ν kinematic viscosity

φ species molar concentration

ρ fluid density

τ relaxation time

0 reference value

1 core region radius

2 annular region inner radius

A annular region

av average quantity

C core region

eq equilibrium conditions

F flow model

(g) gas phase

i coordinate index

j coordinate index

M mass transfer model

n inward normal to solid boundary

r radial direction

(s) solid phase

w value on solid boundary (wall)

z axial direction

1 Introduction

The simultaneous convective-diffusive transport of chemical species and their precipitation onto surfaces following reaction is relevant to many naturally-occurring and engineered systems. Significant amounts of precipitation can lead to complex changes in the boundary geometry, which in turn affects the subsequent fluid flow and species transport. The development of accurate and computationally-efficient models for simulating flow and transport in reacting systems with complex and time-varying domain boundaries is therefore valuable in many applications.

Lattice Boltzmann methods (LBM) have been developed for fluid flow [He and Luo(1997)], heat and mass transport [He, Chen and Doolen (1998)], as well as other physical phenomena such as acoustics and ion transport [Li and Shan (2011); He and Li (2011)].Due to the local nature of the algorithm, lattice Boltzmann methods provide significant advantages in simulating systems with complex and changing geometries, when compared to the finite-difference solution of the Navier-Stokes (N-S) equations, which involves the discretization of numerous derivatives in space. Several LB models have been demonstrated which use one set of distribution functions to determine the velocity field,and an additional set to solve the convective-diffusion equation (CDE), which is coupled to the velocity field. These models can both be implemented with the single relaxation-time

Bhatnagar-Gross- Krook (BGK) collision operator [He, Chen, and Doolen (1998);Parmigiani, Huber, Chopard et al. (2009); Chai and Zhao (2013)], as well as with multiple relaxation times [Yoshida and Nagaoka (2010)].

The LBM for the CDE has been applied to many problems wherein it is desirable to exploit its benefits regarding irregular and/or moving boundaries. Kang et al. simulated the growth of crystals using the LBM to solve the N-S equations and the CDE with a first-order kinetic boundary condition [Kang, Zhang, Lichtner et al. (2004)]. The boundary condition was implemented using lattice-sized control volumes at the liquid/solid interface. Similar approaches were later used to model snow crystal growth in clouds [Lu, Depaolo, Kang et al. (2009)] as well as the growth of hydrate crystals in geological CO2sequestration[Kang, Lichtner, Viswanathan et al. (2010)], and surface growth in reactive capillarydriven flow [Sergi, Grossi, Leidi et al. (2014)]. A similar approach was also used to model flow-related clotting in an investigation of blood clots [Harrison, Smith, Bernsdorf et al.(2007)], however, the passive scalar was treated with a first-order upwind scheme rather than the LBM. These studies used a “stair-case” approximation of the solid boundary, i.e.one in which the boundary does not cut through adjacent lattice boundaries, but rather is staggered between regular lattice nodes, resulting in a “pixelized” boundary. Some other studies have employed a sub-grid representation of the solid boundary, i.e. one not conforming to lattice boundaries. A lattice Boltzmann approach with an immersed subgrid boundary was recently used to simulate solid/liquid phase change [Huang and Wu(2014)]. Here, the curved boundary between the solid and liquid phases was approximated as piecewise linear between a set of Lagrangian points. In another investigation, the growth of dendrite formations in channel flow was simulated using a lattice Boltzmann model for the flow and a phase-field method for the combined mass transfer and solid boundary growth [Hawkins, Angeluta, Hammer et al. (2013)].

Although many applications of these methods occur in cylindrical geometries, i.e. flow in pipes and capillaries, much of the previous work on precipitation/dissolution models has assumed a Cartesian coordinate system. To the authors’ knowledge, there has been no examination to date of problems involving precipitation in cylindrical geometries despite their scientific and industrial relevance. In the current work, we apply recently advanced multiple-relaxation-time lattice Boltzmann models for both incompressible fluid flow and mass transport in axisymmetric cylindrical coordinates. We also employ a new treatment for the third-kind boundary condition for mass transfer on curved boundaries.

The goal of our current work is to develop a basis for predicting the time- and axiallyevolving profile of solid ZnO deposits in a non-isothermal tubular reactor designed for hydrogen production by the heterogeneous hydrolysis of steam with Zn vapor.Heterogeneous hydrolysis with Zn vapor offers a method of water-splitting with higher theoretical efficiency and reliability than previous aerosol-based reactors for hydrolysis with Zn in the Zn/ZnO solar thermochemical cycle [Lindemer, Advani and Prasad (2016);Lindemer, Advani and Prasad (2017)]. The precipitation of solid ZnO during the hydrolysis reaction presents a unique consideration for the design and modeling of reactors for this process. Thus, characterizing the effect of flow, mass transfer, and reaction conditions on the transient accumulation of solid ZnO deposits are the primary reasons for developing the model. In addition, the numerical methods and results presented in this paper may also be useful in understanding reactive precipitation/dissolution processes in other industrial, environmental, and biological flows.The paper is organized as follows: In Section 2, we describe the mathematical models for flow and mass transport that are used in this work. Section 3 presents a detailed description of the assumptions and boundary conditions used for our simulations. Section 4 presents several validation studies that were used to verify the accuracy of the current LB models as well as the associated boundary condition implementations. In Section 5, we present results and discussion from a parametric study to investigate the effects of the relevant non-dimensional parameters (Reynolds, Schmidt, and Damk?hler numbers), other modelspecific parameters, and axially-varying kinetics on the evolution of the ZnO precipitation profile. Finally, in Section 6, we present our conclusions.

2 LB models for flow and mass transport

2.1 Multiple-relaxation-time (MRT) model for axisymmetric incompressible fluid flow

The D2Q9 scheme for axisymmetric incompressible fluid flow presented in Zhou [Zhou(2011)] is used in the current work. This scheme replicates the incompressible Navier-Stokes equations in axisymmetric cylindrical coordinates, which can be written in indicial notation as:

With

by substituting the axisymmetric continuity equation:

into the momentum equation. Here, the momentum equation has been written in a “pseudo-Cartesian” form with the source term SFcontributing the remaining terms in the axisymmetric momentum equation.

Here, the uirepresent the r (radial) and z (axial) components of the hydrodynamic velocity,ρ is the fluid density, P is the pressure, ν is the kinematic viscosity, and repeated indices imply summation over r,z. Li et al. [Li, Mei and Klausner (2013)] extended the BGK scheme presented in Zhou [Zhou (2011)] for implementation with the commonly used MRT collision operator presented in Lallemand et al. [Lallemand and Luo (2000)]. Using this approach, the collision process is described by:

where |f is the vector of pre-collision distribution functions (f0, f1, ... , f8),is similarly the post-collision vector, |feqis the equilibrium distribution vector,

ΩF=diag(1/τF0, 1/τF1, ... , 1/τF8) is a matrix of relaxation parameters, and MFis the collision matrix given by:

The equilibrium distribution functions feqare defined by:

where ρ0is the average density, taken to be 1.0, and c is the sound speed, also taken to be 1.0. The components of the source term vector G are given by:

which correspond the source term SFshown in Eq. (1), as shown by the Chapman-Enskog analysis in Zhou [Zhou (2011)].

After the collision step, the post-collision distributions are then streamed to neighboring nodes according to:

where the (r, z)-components of the discrete velocities eαare defined by:

The weight functions wαare given by w0=4/9, wα=1/9 for α=1-4 and wα=1/36 for α=5-8.Additional details, including the derivation of the scheme for the Cartesian BGK case, can be found in He et al. [He and Luo (1997)]. The kinematic viscosity ν is given in terms of the hydrodynamic relaxation time and the sound speed c as:

where τF=τF7=τF8to enforce isotropic momentum diffusion [Lallemand and Luo (2000)].

After evolving the mesoscopic distribution functions throughout the domain, the local density is then given in terms of the distribution functions by:

the local pressure P is given by:

and the velocities uiare given as:

The hydrodynamic quantities are then used to find the equilibrium distributions fαeqat the next time step.

2.2 Multiple-relaxation-time (MRT) model for mass transport

The D2Q5 model for the passive scalar presented in Li et al. [Li, Mei and Klausner(2013)] is implemented in the current work. The passive scalar in this case is the molar concentration of the reactive species (Zn vapor). The concentration of Zn vapor is assumed to be sufficiently dilute such that the passive scalar approach is accurate.

This scheme replicates the axisymmetric CDE given by:

with

where (ur, uz) are the hydrodynamic velocity components determined by the D2Q9 model presented in Section 2.1, and D is the diffusion coefficient of the reactive species,assumed to be isotropic. Similar to the axisymmetric momentum equation, the CDE has been written in “pseudo-Cartesian” form with the source term SMcontributing the axisymmetric terms.

The collision step for the mass transport model is described by:

where g is the pre-collision vector (g0, g1, ... , g4), g? is similarly the post-collision vector,geqis the equilibrium distribution vector, ΩM= diag(1/τM0, 1/τM1, ... , 1/τM4), and MMis the collision matrix given by:

and the (r, z)-components of the discrete velocities are given by:

It should be noted that the collision matrix MMhas been re-arranged from the definition given in Yoshida et al. [Yoshida and Nagaoka (2010)] and [Li, Mei and Klausner (2013)]for consistency with the numbering of the discrete velocities eαused in the current work. The off-diagonal components of the relaxation matrix are assumed to be zero because diffusion is considered to be isotropic.

The components of the source term vector R are given by:

The equilibrium distributions are given by:

The post-collision distributions then stream to neighboring nodes at the next time step, as in the flow model, according to:

The concentration is given by:

The diffusivity is related to the species relaxation time, time step, and sound speed by:

where τM= τM2= τM4to enforce isotropic diffusion [Zhou (2011)].

3 Model description

3.1 Flow boundary conditions

We consider a core-annular flow configuration as shown in Fig. 1 in which Zn vapor is fed to the central core and steam is fed to the annular region. At the inlet, the velocity boundary condition (given in lattice units) in the core region (0

Similarly we assume fully-developed flow through an annular duct as the inlet boundary condition for R2< r < R:

Figure 1: Schematic of tube reactor. The wall of the Zn-vapor supply tube shown in grey separates the core and annular flow regions, and the modeled region is enclosed by the blue rectangle

where uCand uAare the characteristic velocities assigned in the core and annular regions,respectively. At the reactor inlet, the axial velocity uzis assumed to be zero in the region R1

At the outlet (z=L), a zero gradient outflow condition is assumed for both axial and radial velocity:

The velocities at both the inlet and outlet are specified using the method of Zou et al.[Zhou and He (1997)].

Along the z-axis (r=0), a symmetry condition for the flow velocities is enforced by setting f2=f4, f6= f7, and f5= f8, as discussed in Succi [Succi (2009)]. Thus, the bottom row of lattice nodes resides on the line r=0. All terms with a factor of 1/r in the collision step are neglected for these nodes, as they become zero when L’H?pital’s rule is applied at r=0 [Premnath and Abraham (2005); Reis and Phillips (2008); Zhou (2011)].

The flow boundary condition on the solid boundary is implemented according to the interpolated bounce-back method for curved/moving boundaries presented in Lallemand et al. [Lallemand and Luo (2002)].

Along the tube wall, the growth of the solid precipitate is assumed to progress in the direction normal to the existing surface. The hydrodynamic velocity boundary condition for the solid boundary is determined by the growth rate of the solid layer projected onto each coordinate direction, similar to the condition in Li et al. [Li, Huang and Meakin(2010)]:

where φ(s)is the constant molar concentration of the bulk solid, which is taken to be 1000.0 unless otherwise stated, r’’ is the local reaction rate, and θw(z) is the local angle of the curved solid boundary measured from the z-direction. The local angle of the wall is determined from the current shape of the solid deposit which is described by the curve r"=W(z),

When enforcing the bounce-back condition, the value of θwat the point of the outgoing population’s intersection with the solid boundary is determined using quadratic Lagrangian interpolation at nearest-neighbor values of z, (i.e. at z-1,z,z +1), as the intersection point will not generally be on a node. The distance between a boundary node and the boundary (i.e. between xfand xffin Fig. 2) is determined by using a piecewiselinear representation of the boundary curve r=W (z), so that the computation is simple.This piecewise linear description is similar to that used in Huang et al. [Huang and Wu(2014)], however rather than using a dynamically-updated set of Lagrangian points to represent the interface, here, the curve is tracked using a simple Eulerian description, i.e.the curve’s height is tracked at a fixed set of z-values.

Figure 2: Schematic of boundary treatment at solid-fluid interface

3.2 Mass transfer boundary conditions

The inlet concentration boundary condition in the core region (0≤r≤R1) is:

and φ = 0.0 for r > R1. The value of φ0is chosen so that the baseline molar concentration of Zn vapor appropriately lower than the chosen molar concentration of the solid phase.The Dirichlet boundary condition for concentration is implemented using the “scheme D”presented in Liu et al. [Liu, Lin, Mai et al. (2010)]. At the outlet, a zero-gradient condition is assumed for the concentration:

which is implemented using the method presented in Liu et al. [Liu, Lin, Mai et al. (2010)]for Neumann boundary conditions. Similar to the flow model, a symmetry boundary condition is implemented on r=0 by setting g2=g4. The collision terms involving 1/r are likewise ignored on these nodes.

The reaction rate is assumed to follow a first-order kinetic model in terms of φw, the concentration on the wall:

We also assume that the mole balance in the normal direction on r=W (z) includes the mass flux induced by the growth of the solid layer, as in Li et al. [Li, Huang and Meakin(2010)]:

where n is the normal direction pointing into the fluid domain, and k is the kinetic constant. In order to implement this boundary condition on a curved boundary, we used the method in [Li, Mei and Klausner (2013)], in which both the tangential and normal mass fluxes are calculated at the intersection of the outgoing population eαand the boundary, and are then projected onto the incoming population with velocity eαˉ. The presence of a tangential mass flux on a curved boundary is an important consideration for accurate implementation of this boundary condition.

Fig. 2 illustrates how the boundary condition in Eq. (29) is implemented. The values of g?αand g?αˉat points xfand xffare used in the calculation as well as g?βˉ and g?βat xf′and xf′f. The vectors eαand eβare perpendicular, and constitute basis vectors for calculation of the tangential and normal mass fluxes. The green points are used to interpolate the value of g?β(xf′f), which generally does not reside on a regular node point. Values at xf′are determined from extrapolation from xfand xff, as described in Guo et al. [Guo, Zheng,and Shi (2002)]. The curved boundary, shown in blue, is approximated as piecewise linear,as shown by the red lines.

The normal and tangential mass fluxes are projected onto the αˉ-direction using the known relationships between θw, θnαˉ, and θnβˉ. Finally, the projected mass flux is implemented as a Neumann boundary condition, as described in Li et al. [Li, Mei, and Klausner (2013)].

Once all incoming populations have been determined, the wall concentration profile φw(z) is then calculated from g2on the boundary nodes using Eq. (41b) in Li et al. [Li,Mei and Klausner (2013)]. A profile of the reaction rate can then be determined from the definition r′′(z) = kφw(z).

3.3 Non-dimensional parameters

The relevant non-dimensional numbers in this system are the Reynolds number, given by:

with

where QAand QCare the volumetric flow rates through the annular and core regions of the inlet, respectively, determined analytically from Eqs. (20) and (21). Nris the number of lattice units in the radial direction.

The Mach number is defined as:

The Schmidt number is defined as:

The Peclet number is defined as:

The Damk?hler number is defined as:

We also define a second Damk?hler number as the ratio of the reaction rate to the advective mass transfer rate:

In the current work, this parameter is only of interest for cases with axially-varying kinetics, where k1represents the initial value of the kinetics constant at the inlet. The non- dimensional time is defined using diffusive scaling:

where Ntis the number of time steps.

4 Model validations

To verify the accuracy of the models used, several validation problems were selected. First,to validate the coupling of the flow model to the convective-diffusion equation, the problem of hydrodynamically-developed but thermally-developing flow, also known as the Graetz problem, was simulated using the LBM code. At the outset, a uniform body force was applied to the flow model in a circular tube with periodic boundary conditions at the inlet and outlet, and the flow was allowed to converge to a fully-developed Poiseuille velocity profile. Next, a constant temperature of T0=1.0 was applied at the inlet and a constant temperature of Tw=0.0 was applied at r=R, with a symmetry boundary condition at r=0.The Nusselt number is defined as Nu=(R), and θ=T/T the bulk temperature0θ=∫R2πruT dr /T∫R2πrudr. The temperature gradient at the wall was determined byb000 non-equilibrium extrapolation and the bulk temperature was obtained by numerical integration of the temperature and velocity profiles. The results for Re=250 and Pr=0.7 are shown in Fig. 3 along with the theoretical results from the first 10 terms of the Graetz series solution. The agreement with theory is found to be quite good except for very low values of z due to the singular point at (z=0,r =R). The Nusselt number closely approaches the well-known result of 3.657 at the outlet. In addition, the axial evolution of the bulk temperature computed by LBM matches very well with theory.

Figure 3: Results for thermally-developing flow with a constant wall temperature boundary condition for Re = 250, Pr = 0.7

Second, to validate the movement of the solid boundary as well as the implementation of the flow and mass transfer boundary conditions, a one-dimensional reaction-diffusion validation problem similar to the one given in Li et al. [Li, Huang and Meakin (2010)]was simulated. In this problem, a solute diffuses and reacts on a wall, which moves in accordance with the growth of the deposited solid over time. This problem is formulated in Cartesian coordinates, so all of the axisymmetric terms were left out of the collision operator for mass transfer in Eq. (14). The problem is implemented with five lattice points in the wall-parallel (x) direction, and periodic boundary conditions in x so that the problem has no x-dependence. Initially, the entire domain has a uniform concentration φ0; subsequently, the concentration at the y-origin (y=0) remains fixed at φ0, so that φ(0, t) = φ(y, 0) = φ0. The initial location of the wall is given by W(t=0)=W0. Assuming a first-order reaction rate, the velocity of the wall in the wall-normal direction is proportional to the reaction rate on the surface (y = W), and is given by:

where φ(W) is the concentration at the wall, k is the kinetic constant, and φ(s)is the constant concentration of the solid phase. Then, for reaction-limited slow growth,(Da<<1, Pe<<1), the transient solution for the wall location is given as:

The LBM results are compared with the analytical solution in Fig. 4 for Ny= 40, with time non-dimensionalized using the diffusive scaling as in Eq. (37). The two are found to diverge slightly at lower values of Was fewer lattice points are being used in the simulation, and hence there is a loss of accuracy with increasing time.

Figure 4: Comparison of LBM results with theory for the one-dimensional reactiondiffusion problem under slow growth conditions

Finally, as a check on the transient performance of the axisymmetric model and overall convergence, the problem of transient diffusion in slug flow with sinusoidally-varying wall concentration, as described in Li et al. [Li, Mei and Klausner (2013)] was simulated. We consider a section of a pipe of radius R and length L=2R with a uniform axial velocity U. The concentration boundary condition is prescribed as φ(R,z,t)=cos(kz + ωt) with k=2π/L and ω=2π/Γ, where Γ is the time period of oscillation. The Strouhal and Peclet numbers are defined as St=and Pe=respectively. Periodic boundary conditions are applied at z=0 and z=L.

The analytical solution for the concentration in the pipe is given by Li et al. [Li, Mei and Klausner (2013)]:

The convergence of the temporally and spatially averaged error E2is shown in Fig. 5 for different values of the mass transfer relaxation time τM, for Pe=20.0, St=10.0. The convergence is seen to be second-order with the number of lattice points in the radial direction (Nr) as was also found in Li et al. [Li, Mei and Klausner (2013)].

Figure 5: Normalized cumulative global error E2 after one period Γ for St=10, Pe=20

5 Results and discussion

Following the validation studies presented in the previous section, the model is now applied to investigate the heterogeneous hydrolysis of steam with Zn vapor in a tubular reactor under a negative axial temperature gradient. The negative axial temperature gradient is of interest for optimization of the hydrolysis reaction in the Zn/ZnO thermochemical cycle for water-splitting [Lindemer, Advani and Prasad (2016)]. For all cases, the simulation is initialized with the fluid at rest and the Zn vapor concentration set to zero everywhere. Thus,initially, only the background flow of steam is present. The flow boundary conditions are then applied and the flow is allowed to develop until it reaches a steady-state, and then the mass transfer boundary conditions are applied. This sequence is followed in order to decouple the development of the initial flow field from the results for different cases. The simulations use Nr=40 for all cases, which was found to give sufficient accuracy with reasonable computational time in the model validations.

The evolution of the axial velocity uzand Zn-vapor concentration φ for a typical case with Sc=0.7, Da=1.0, and Re=1.0 is shown in Figs. 6 and 7, respectively for three instants in time. As the solid ZnO layer (shown by the dashed white line) grows inward, the axial velocity in the throat region increases since the inlet mass flux is fixed. The velocity boundary layer is also seen to be thinnest at the point of fastest deposition since the reduction in cross-sectional area causes the flow to accelerate. By numerically integrating the velocity profile along the radial direction, it was verified that the mass flow rate remains constant for all z locations, with any small differences owing to different numbers of radial nodes at each value of z resulting in different degrees of error.

Figure 6: Contours of uz (lattice units) for Sc=0.7, Da=1.0, and Re=1.0, at t?=612.5,6862.5, 9362.5. The dashed white line indicates the deposition profile at each time step

Figure 7: Contours of φ for Sc=0.7, Da=1.0, and Re=1.0, at t?=612.5, 6862.5, 9362.5.The dashed white line indicates the deposition profile at each time step

Due to the relatively low value of Da in this case, the concentration contours in Fig. 7 do not show the presence of a boundary layer adjacent to the solid ZnO boundary. It is also apparent that the contours of constant concentration are stretched axially at later time steps owing to the stronger effect of advection due to increasing velocities. For increasing values of Da, the results are qualitatively similar to Fig. 7, except with a more pronounced concentration boundary layer, as the case of Da>>1 represents a zero-concentration boundary condition.

5.1 Effect of Re

Fig. 8 shows the effect of increasing Re for Re 1.0 for Da=20.0, Sc=0.7, and Qr=1.0 for different time steps. The plots for Re=0.01 and Re=0.1 are virtually identical since for very low flow velocities the problem is essentially diffusion-dominated with advection playing a minimal role in mass transport. For Re=1.0, the deposit exhibits a similar profile to the lower Re cases, however the growth is faster over each time interval. This is because the φu term in the total inlet mass flux is larger at higher flow velocities, and hence the total mass flux into the system is significantly larger for Re = 1.0 compared to the lower Re cases.

Figure 8: Deposition profiles for different values of Re, with Da=20.0, Sc=0.7, QA/QC=5.0. Results for t?=508.0 are shown in blue, t?=848.8 in green, and t?=1530.7 in red

Fig. 9 shows the effect of Re on the evolution of the deposition profile for Re>1, with all other non-dimensional quantities the same as those in Fig. 8. Increasing Re primarily has the effect of changing the shape of the deposition profile such that the axial location of fastest deposition is shifted downstream, as well as increasing the total amount of deposition over any time interval. The shifting of the deposition profile downstream with increasing Re is due to the advection of mass increasing downstream with Re, and the increased rate of overall deposition is again due to higher advective mass flux into the system at higher Re. Hence, for constant Sc, the value of Re is a strong predictive factor in the problem, with Re<1 representing one regime with minimal differences in the evolution of the ZnO deposition profile, and Re>1 representing a regime with significant sensitivity to Re.

Figure 9: Deposition profiles for different values of Re, with Da=20.0, Sc=0.7, QA/QC=5.0.Results for t?=167.0 are shown in blue, t?=508.0 in green, and t?=848.8 in red

5.2 Effect of Da

Fig. 10 shows the effect of Da for moderate-to-high values of Da with Re=10.0, Sc=0.7,and Qr=5.0. At low t?, the deposition profiles for Da=5.0, 10.0, and 20.0 are nearly identical. This is because for this range of Da, the kinetics are very fast relative to the rate of diffusion, and hence the reaction is mass transport-limited. At higher values of t?, the radial diffusion length decreases due to the growth of the solid deposit, and thus the mass transfer resistance due to diffusion also decreases. In this regime, the kinetics begin to play a stronger role in the rate of reaction, and the deposition profiles become more sensitive to Da, with higher Da resulting in faster deposition. Faster deposition results in a faster decrease in diffusion length, and thus the differences in the deposition profiles for different Da become more pronounced over time. At high Da the axial location of fastest deposition is rather insensitive to Da, because the initial location is strongly mass transfer-controlled.

Fig. 11 shows the effect of Da for Da=1.0, with all other conditions the same as those in Fig. 10. In this range, the reaction is kinetics-limited, and hence the deposition profiles immediately show sensitivity to Da at low t?. For higher Da, the deposition is faster and occurs further upstream, as a higher amount of reactive species that contacts the wall will react rather than be transported downstream. The long vertical tick marks in Fig. 11 indicate the axial position of fastest deposition for each case. Thus, while increasing Re has the effect of pushing the location of fastest deposition downstream for constant Da,increasing Da has the opposite effect for constant Re.

Figure 10: Deposition profiles for different values of Da, Re=10.0, Sc=0.7, QA/QC=5.0

Figure 11: Deposition profiles for different values of Da, Re=10.0, Sc=0.7, QA/QC=5.0.Results for t?=249.5 are shown in blue, t?=758.8 in green, and t?=1268.1 in red

5.3 Effect of Sc

In Section 5.1, the effect of changing Re for constant Sc was examined, which is equivalent to changing Pe. The effect of changing Sc and Re simultaneously for Pe=10.0 and Da=20.0 was investigated. It was found that when Pe is held constant, the results are virtually identical for any value of Sc. Any slight differences in the deposition profiles can be attributed to differences in numerical error, as the value of M a is different for each case. This insensitivity to Sc for constant Pe indicates that in the moderate Re regime, the actual predictor of the deposition profile is Pe. However, it is possible that under turbulent flow conditions, Re would play a role independent of the value of Pe as discussed in Hawkins et al. [Hawkins, Angeluta, Hammer et al. (2013)], but an investigation into this effect is beyond the scope of the current work.

5.4 Effect of QA/QC

The effect of changing the annular-to-core flow rate ratio QA/QCwas also investigated while holding all other non-dimensional parameters constant. Despite the non-uniform velocity profiles at the inlet, for a constant value of Re, the hydrodynamic entry length is roughly the same for each case, and hence the axial location of fastest deposition is nearly identical for all cases examined. This is likewise found to hold for lower values of Da.Although the inlet velocity profiles are very different for these cases, the flow develops to nearly the same velocity profile at a relatively low value of z/R for all cases. As the ratio QA/QCis increased for constant Re, QCis decreased, which decreases the total flux of reactive species into the system. For higher QA/QC, the total amount of deposition is decreased for any value of t?. Thus, the value of QCis important in determining the total mass flux of reactant into the system, and thus the overall rate of deposition, but the results are relatively insensitive to QA/QCfor any given QC.

5.5 Effect of φ(s)

Here, we examine the effect of the product molar density φ(s) on the deposition profile.As discussed in Li et al. [Li, Huang, and Meakin (2010)], the interface velocity uwis typically very small for gas/solid precipitation, and in fact can be safely ignored for large values of φ(s) with minimal impact on the results. Thus, for high values of φ(s), the problem is quasi-steady in time. Deposition profiles are shown for different values of φ(s)in Fig. 12. Physically, different values of φ(s) represent different solid products, or different morphologies of the same material. Hence, a product with a higher porosity would correspond to a solid with a lower molar density. The deposition profile for φ(s)=500.0 at the earliest time overlaps with the profile for φ(s)=1000.0 at the second time instant, and with the profile for φ(s) =2000.0 at the fourth time instant. This is to be expected, as the growth of the deposit should be twice as slow when the density of the solid is doubled, and so on. This illustrates that the shape of deposition profile for a material with a given φ(s) can be easily inferred from results for other values of φ(s), as the results scale predictably for slow growth conditions. As φ(s) is decreased significantly for a given k and φ0, the wall velocities will become comparable to the axial velocity of the flow, and the growth may not be quasi-steady; however, investigation of this regime this is beyond the scope of the current work.

Figure 12: Deposition profiles for different values of φ(s), Re=10.0, Sc=0.7, Da=10.0, and QA/QC=5.0

5.6 Effect of Da2 with variable kinetics

The heterogeneous hydrolysis of steam by Zn vapor is optimally performed nonisothermally under a negative axial temperature gradient, as discussed in (15) and experimentally verified in Lindemer et al. [Lindemer, Advani, and Prasad (2017)]. The reaction should ideally begin at a high temperature to minimize the use of carrier gas, and the reactor temperature should then decline downstream in order to take advantage of faster kinetics and obtain higher equilibrium yields. Thus, it is of particular interest to examine cases where the kinetic constant varies due to a negative temperature gradient in the axial direction.

Accordingly, we now examine the effect of Da2=k1/uavfor cases with axially-varying kinetics. The kinetic constant is assumed to be constant at k=k1for zz0. For the following cases, we also define Da=k1R/D in order to characterize the initial reaction rate separately from the axial variation in kinetics. The kinetic constant is used as a proxy for the tube wall temperature, and we assume that the temperature gradient does not affect the flow or mass transfer. The increase in k is used to simulate both the increase in kinetics and the shifting of equilibrium conditions toward ZnO production at lower temperatures [Lindemer, Advani and Prasad (2016)],although the current rate law does not explicitly account for reaction equilibrium.

Results are shown in Fig. 13 for Da=0.1, k2/k1=50.0, Sc=1.0, z0=5Nr, and various values of Re. Since Da2can be written as Da2=Da/ReSc, we vary Re in order to change Da2while holding the kinetics profile, Sc and Da constant. As Re is increases, the deposition profile transitions from exhibiting a single maximum to exhibiting two local maxima. This phenomenon was also consistently observed in experimental results that were obtained using a non-isothermal laboratory-scale tube reactor for the heterogeneous hydrolysis of steam with Zn vapor [Lindemer, Advani, and Prasad (2017)]. Moreover, as Re increases,the relative size of the second maximum increases and the location of the maximum occurs further downstream.

Figure 13: Deposition profiles for the varying kinetics case with Sc=1.0, Da=0.1, QA/QC=5.0, and various values of Re. Results for t?=1516.6 are shown in blue, t?=7707.1 in green, and t?=13897.6 in red

For a given kinetics profile, the appearance of this two-peak deposition profile is found to depend on Da2, which is the ratio of the reaction rate to the advective transport rate in the constant kinetics (entrance) region. As Re is increased with constant Da and Sc, Da2is decreased by definition, and the relative size of the second deposition peak is found to increase, as shown in Fig. 13. This is because as the flow rate increases, the kinetics at the inlet become slower compared to advection and thus more species is advected downstream rather than deposited. The concentrations are thus higher in the downstream section with faster kinetics, resulting in a rapid increase in reaction rates. For high enough values of Da2, the second peak in the deposition profile does not occur, as very little reactant remains to be transported to the downstream region with faster kinetics. It is also found that for a constant value of Da2, varying the value of Sc does not significantly change the shape of the deposition profile, and therefore has no effect independent of Da2, similar to what was found for the constant kinetics case in Section 5.3.

For a constant value of Da2, the shape of the kinetics profile also influences the shape of the deposition profile. Figs. 14 and 15 show the effect of varying the values of z0and k2/k1,respectively for Da2=0.1. As shown in Fig. 14, as z0is increased, the second increase in the deposition is pushed correspondingly further downstream, shortly downstream from z0.Increasing the value of k2/k1increases the relative size of the second maximum, as shown in Fig. 15. At a given Da2and z0, for low enough values of k2/k1, the second maximum will not occur, as the slow increase in kinetics does not make the case significantly different from a constant kinetics case, as shown in Fig. 15 for k2=20k1. As z0→0, only a single deposition maximum is observed, because most of the limiting reactant is depleted very close to the entrance due to the sharply increasing kinetics. Likewise, if z0is too large, only one maximum is observed, because most of the limiting reactant is depleted before z0is reached.As Da2is increased, it is found that a larger value of k2/k1is required for the same value of z0, in order for the two-maxima deposition profile to occur. This is because the effect of advection is weaker in the constant kinetics section, resulting in lower concentrations in the high kinetics region, which must then be overcome with a sharper increase in kinetics.Similarly, as Da2is increased, a lower value of z0is required to produce a two-maxima deposition profile for the same value of k2/k1, as the kinetics must rise sharply before a significant amount of the limiting reactant has been depleted.

Figure 14: Deposition profiles for the varying kinetics case with Re=1.0, Sc=1.0, Da=0.1,QA/QC=5.0, and various values of z0. Results for t?=1516.6 are shown in blue, t?=7707.1 in green, and t?=13897.6 in red

6 Conclusions

The reactive gas/solid precipitation of ZnO in a circular tube has been simulated using multiple-relaxation-time lattice Boltzmann schemes in axisymmetric cylindrical coordinates for both flow and mass transfer. A D2Q9 scheme is used for the incompressible fluid flow, and a D2Q5 scheme is used for the mass transfer. The model assumes first-order kinetics at the solid boundary and a third-kind boundary condition for curved surfaces. The flow and mass transfer models are both validated using several problems in axisymmetric cylindrical coordinates with analytical solutions.

The deposition of solid products is found to depend strongly on whether the mass transfer is diffusion or advection-limited, with the results being nearly identical for diffusion limited conditions due to the assumption of constant concentration at the inlet. In the advectionlimited regime, the axial location of fastest mass deposition is found to scale with Pe, or equivalently with Re for constant Sc. The results are insensitive to Sc for a given value of Pe. It is also found that the results depend strongly on whether the reaction is kinetics- or diffusion-limited, with the axial location of fastest deposition scaling inversely with Da under diffusion-limited conditions. Under diffusion-limited conditions, the reaction transitions from being diffusion-limited to kinetics-limited as the solid deposit grows and diffusion lengths are decreased. Thus, over short time scales, the results are nearly independent of Da for Da>5.

Finally, the effect non-isothermal conditions are simulated by axially varying the kinetic constant. The results show that the value of Da2=Da/ReSc determines whether the deposition profile has one or two peaks for a given kinetics profile, with the two-peak pattern occurring at lower values of Da2. Likewise, the shape kinetics profile has an effect on the deposition profile, with sharp increases in kinetics at moderate distances from the tube entrance tending to produce the second maximum in the deposition profile.

The model and results help to gain insight into the reactive deposition of ZnO in heterogeneous water-splitting using Zn vapor. More generally, this approach may be applicable to similar heterogeneous reactive/precipitation processes in axisymmetric cylindrical coordinates in other engineering situations. This work has demonstrated methods of validating such a model, and discusses the key physics involved in this process as elucidated from our numerical results.

Acknowledgement:Funding for this work was provided by the Federal Transit Administration, the University of Delaware’s NSF-IGERT Sustainable Energy from Solar Hydrogen program, and the Air Products and Chemicals, Inc. Ph.D. Fellowship.