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

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 temporallyand 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.


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 CO 2 sequestration [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.

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: 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 S F contributing the remaining terms in the axisymmetric momentum equation.
Here, the u i represent the (radial) and (axial) components of the hydrodynamic velocity, is the fluid density, is the pressure, is the kinematic viscosity, and repeated indices imply summation over , . ] 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 | is the vector of pre-collision distribution functions ( 0 , 1 , ... , 8 ), |^ is similarly the post-collision vector, | is the equilibrium distribution vector, F =diag(1/τ F0 , 1/τ F1 , ... , 1/τ F8 ) is a matrix of relaxation parameters, and F is the collision matrix given by: The equilibrium distribution functions are defined by: where ρ 0 is 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 are given by: which correspond the source term S F shown 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: The weight functions w α are given by w 0 =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 =τ F8 to 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: (11) the local pressure P is given by: 12) and the velocities u i are given as: The hydrodynamic quantities are then used to find the equilibrium distributions at the next time step.

Multiple-relaxation-time (MRT) model for mass transport
The D2Q5 model for the passive scalar presented in Li et al. ] 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: where (u r , u z ) 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 S M contributing the axisymmetric terms. The collision step for the mass transport model is described by: where g is the pre-collision vector (g 0 , g 1 , ... , g 4 ), ĝ is similarly the post-collision vector, g eq is the equilibrium distribution vector, M = diag(1/τ M0 , 1/τ M1 , ... , 1/τ M4 ), and is 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 has been re-arranged from the definition given in Yoshida et al. [Yoshida and Nagaoka (2010)] and ] for consistency with the numbering of the discrete velocities α 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 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: (22) The diffusivity is related to the species relaxation time, time step, and sound speed by: [Zhou (2011)].

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<r<R 1 ) is assumed to be fully-developed flow through a circular duct: Similarly we assume fully-developed flow through an annular duct as the inlet boundary condition for R 2 < 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 u C and u A are the characteristic velocities assigned in the core and annular regions, respectively. At the reactor inlet, the axial velocity u z is assumed to be zero in the region R 1 <r<R 2 which corresponds to the wall thickness of the tube that introduces the core flow into the reactor. The radial velocity u r is also assumed to be zero for all r at z=0. 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 2 = 4 , 6 = 7 , and 5 = 8 , 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/ 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, '' 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 " = ( ), When enforcing the bounce-back condition, the value of θ w at the point of the outgoing population's intersection with the solid boundary is determined using quadratic Lagrangian interpolation at nearest-neighbor values of , (i.e. at − 1, , + 1), as the intersection point will not generally be on a node. The distance between a boundary node and the boundary (i.e. between and in Fig. 2) is determined by using a piecewiselinear representation of the boundary curve = ( ), 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 -values.

Mass transfer boundary conditions
The inlet concentration boundary condition in the core region (0≤r≤R 1 ) is: and φ = 0.0 for r > R 1 . The value of φ 0 is 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 g 2 =g 4 . The collision terms involving 1/ are likewise ignored on these nodes. The reaction rate is assumed to follow a first-order kinetic model in terms of , the concentration on the wall: r ″ = kϕ w (31) We also assume that the mole balance in the normal direction on = ( ) includes the mass flux induced by the growth of the solid layer, as in Li et al. [Li, Huang and Meakin (2010)]: where is the normal direction pointing into the fluid domain, and is the kinetic constant. In order to implement this boundary condition on a curved boundary, we used the method in ], 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 and are used in the calculation as well as gˆβ¯ and gˆβ at ′ and ′ . 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ˆβ( ′ ), which generally does not reside on a regular node point. Values at ′ are determined from extrapolation from and , 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 ]. Once all incoming populations have been determined, the wall concentration profile ( ) is then calculated from g 2 on the boundary nodes using Eq. (41b) in Li et al. ]. A profile of the reaction rate can then be determined from the definition ′′( ) = ( ).

Non-dimensional parameters
The relevant non-dimensional numbers in this system are the Reynolds number, given by: where Q A and Q C are the volumetric flow rates through the annular and core regions of the inlet, respectively, determined analytically from Eqs. (20) and (21). N r is 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 k 1 represents the initial value of the kinetics constant at the inlet. The non-dimensional time is defined using diffusive scaling: where N t is the number of time steps.

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 T 0 =1.0 was applied at the inlet and a constant temperature of T w =0.0 was applied at r=R, with a symmetry boundary condition at r=0. The Nusselt number is defined as = −2 ( ), and = / 0 the bulk temperature The temperature gradient at the wall was determined by 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 due to the singular point at ( = 0, = ). 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. 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 ( = 0) = 0 . 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 = ), and is given by: where φ( ) 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 N y = 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.
where = �1 + + 2 , and I 0 is the modified first-kind Bessel function of order 0. The convergence of the temporally and spatially averaged error E 2 is 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 (N r ) as was also found in Li et al. ].

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

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 N r =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 u z and 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.   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. Fig. 8 shows the effect of increasing Re for Re 1.0 for Da=20.0, Sc=0.7, and Q r =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.  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. Results for t * =167.0 are shown in blue, t * =508.0 in green, and t * =848.8 in red Fig. 10 shows the effect of Da for moderate-to-high values of Da with Re=10.0, Sc=0.7, and Q r =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. Results for t * =249.5 are shown in blue, t * =758.8 in green, and t * =1268.1 in red

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.

Effect of Q A /Q C
The effect of changing the annular-to-core flow rate ratio Q A /Q C was 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 Q A /Q C is increased for constant Re, Q C is decreased, which decreases the total flux of reactive species into the system. For higher Q A /Q C , the total amount of deposition is decreased for any value of t * . Thus, the value of Q C is 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 Q A /Q C for any given Q C.

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 u w is 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.

Effect of Da 2 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 Da 2 =k 1 /u av for cases with axially-varying kinetics. The kinetic constant is assumed to be constant at k=k 1 for z<z 0 and then increases as = 1 + 2 � − 0 /2 � 2 for z>z 0 . For the following cases, we also define Da= k 1 R/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, k 2 /k 1 =50.0, Sc=1.0, z 0 =5N r , and various values of Re. Since Da 2 can be written as Da 2 =Da/ReSc, we vary Re in order to change Da 2 while 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. For a given kinetics profile, the appearance of this two-peak deposition profile is found to depend on Da 2 , 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, Da 2 is 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 Da 2 , 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 Da 2 , varying the value of Sc does not significantly change the shape of the deposition profile, and therefore has no effect independent of Da 2 , similar to what was found for the constant kinetics case in Section 5.3. For a constant value of Da 2 , 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 z 0 and k 2 /k 1 , respectively for Da 2 =0.1. As shown in Fig. 14, as z 0 is increased, the second increase in the deposition is pushed correspondingly further downstream, shortly downstream from z 0 . Increasing the value of k 2 /k 1 increases the relative size of the second maximum, as shown in Fig. 15. At a given Da 2 and z 0 , for low enough values of k 2 /k 1 , 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 k 2 =20k 1 . As 0 → 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 z 0 is too large, only one maximum is observed, because most of the limiting reactant is depleted before z 0 is reached. As Da 2 is increased, it is found that a larger value of k 2 /k 1 is required for the same value of z 0 , 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 Da 2 is increased, a lower value of z 0 is required to produce a two-maxima deposition profile for the same value of k 2 /k 1 , as the kinetics must rise sharply before a significant amount of the limiting reactant has been depleted. in green, and t * =13897.6 in red CMES, vol.117, no.3, pp.527-553, 2018 Figure 15: Deposition profiles for the varying kinetics case with Re=1.0, Sc=1.0, Da=0.1, Q A /Q C =5.0, and various values of k 2 /k 1 . Results for t * =1516.6 are shown in blue, t * =7707.1 in green, and t * =13897.6 in red

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 Da 2 =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 Da 2 . 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.