of Chemical Engineering NUMERICAL SIMULATION OF THE INCOMPRESSIBLE TURBULENT FLOW OF A BINARY MIXTURE OF AIR-WATER VAPOR : APPLICATIONS IN DRYING PROCESS

The present work deals with numerical simulation of the incompressible turbulent flow of a binary mixture air-water vapor inside channels. Calculations are performed using the RANS (Reynolds Average Navier-Stokes Equations) formulation in addition to the standart k-ε turbulence model. The mathematical model is discretized by a finite difference scheme, being adopted second order accurate expressions for both convection and diffusion terms. The mesh arrangement is collocated and artificial dissipation terms are added to control the odd-even decoupling problem. The numerical scheme is applied to solve the flow of a binary mixture of air-water vapor inside plane channels and sudden expansions. The validation performed indicates that the present method reproduces satisfactorily the literature data for both concentration profiles and Sherwood number. Furthermore, the parametric analysis performed indicates that the drying process (wall mass flux) is very sensitive to the flow parameters investigated, i.e., inlet flow velocity and channel expansion ratio.


INTRODUCTION
The turbulent flow of binary mixtures occurs in practical problems of the chemical and mechanical engineering areas.In the convective drying process, water vapor is removed from moist objects (e.g.wood, food, etc.), being transported by the hot air stream.It is important to mention that convective drying is the most used method in agricultural, food, paper and timber processing industries (Mohan and Talukdar, 2010).In general applications, convective drying of moist objects conjugates heat and mass transfer both inside the solid and in the flowing medium (Vaz Jr et al., 2013).The understanding of the related physics of the problem is important for industrial process control.
The literature works dealing with this topic are increasing, mainly in the numerical analysis area.Younsi et al. (2010) adopted the commercial software ANSYS-CFX to solve numerically the turbulent flow of a binary mixture in a wood dry kiln.The authors compared the numerical results with experimental data aiming at validating the simulation.Kadem et al. (2011) solved numerically the conjugate heat and mass transfer problem in a convective wood dry kiln using the commercial software FEMLAB.The mathematical model adopted by the authors to solve the flow of a binary mixture is the Navier-Stokes equations subjected to laminar flow conditions.The results obtained indicate a correlation between drying rate in the solid and the flow behaviour.In a similar route, Lamnatou et al. (2009) studied numerically the drying process of a moist object (a porous sliced wood) subjected to incompressible laminar flow of a binary mixture inside a channel.The authors analysed the effects of two flow configurations, and the results indicated that using the upstream flow divider device produced a higher mass transfer rate at the surface of the moist object (enhancing the drying rate).
In other interesting practical applications, some literature works solve the flow of a binary mixture to study both the environmental problem of pollutant dispersion and the chemical corrosion of metallic surfaces (Xiong et al., 2014;Lateb et al., 2013;Arellano et al., 2013;Arellano and Rivera, 2014).The aforementioned works adopted the numerical methodology to study practical aspects of the problem, i.e., to estimate the mass flux of the chemical component on the surface, the concentration of the pollutant between the buildings in a city, etc.
The present work deals with the 2D numerical simulation of turbulent flow of a binary mixture (airwater vapor) aiming to study some aspects of the convective drying process.The numerical scheme applied by the authors was originally proposed for solving turbulent Newtonian flows (Zdanski et al. 2004).The present methodology was also applied successfully to solve non-Newtonian polymer melt flows inside channels (Zdanski et al., 2008;Zdanski and Vaz Jr., 2011).Therefore, the present work performs an extension of the author's methodology for solving turbulent flows of binary mixtures, including the solution of the concentration equation.The results obtained indicate that the present method reproduces satisfactorily the literature data for both concentration profiles and Sherwood number.Furthermore, the parametric analysis performed shows that the drying process (mass flux at the channel walls) is dramatically affected by the flow parameters investigated, i.e., inlet flow velocity and channel expansion ratio.

Governing Equations
The mathematical model employed corresponds to the classical Reynolds Averaged Navier-Stokes Equations (RANS) together with the Boussinesq eddy viscosity approximation.Within this framework, the mass, linear momentum, energy and chemical species conservation laws are solved to obtain the average velocities, pressure, temperature and concentration profiles of the flow.In cartesian tensor notation, the equations are given by (Kays and Crawford, 2005): the turbulent eddy viscosity being calculated by: Brazilian Journal of Chemical Engineering Vol. 33, No. 02, pp. 287 -296, April -June, 2016 The constants appearing in the preceding equations are empirically determined (Launder and Spalding, 1974), i.e., 0.09 The turbulent diffusion coeficients in the energy and chemical species equations are obtained by using the analogy between linear momentum, heat and mass diffusion processes through the definition of the turbulent Prandtl and Schmidt numbers (Kays and Crawford, 2005), The turbulent Prandtl and Schmidt numbers are assumed to be constants, being adopted common values for engineering calculations Pr 0.9 and Crawford, 2005).It is worth to mention that the standard    turbulence model is valid only for the fully turbulent region; in the near wall regions the law of the wall for velocity, temperature and concentration profiles is adopted (Launder and Spalding, 1974;Kays and Crawford, 2005;Arpaci and Larsen, 1984).

Numerical Methodology
The mathematical model is discretized with a finite difference scheme, second order accurate formulae being adopted for both the convection and diffusion terms.The mesh arrangement is collocated and artificial dissipation terms are added to control the odd-even decoupling problem.The pressure-velocity coupling strategy employed comprises the solution of the Poisson equation for pressure.The solution procedure follows a pseudo-transient march in time, aiming at the steady-state solution.The present work constitututes an extension of the author's methodology (Zdanski et al., 2004;Zdanski et al., 2008;Zdanski and Vaz Jr., 2011) for solving turbulent flows of binary mixtures.The main general steps of the present scheme are given in sequence.
The linear momentum, energy and chemical species conservation equations, as well as the turbulence model equations, can be written in a compact form (for a 2D cartesian coordinate system) where Q is the vector of conserved variables, E, F and S are the flux vectors containing the convection, diffusion and source terms.The vector Q is expanded in time by Taylor expansion: in which m and t  are the time step and time incre- ment, respectively.Furthermore, from Equation ( 11) one finds: Solution of Equation ( 12) requires linearization of vectors E, F and S in terms of vector Q using Taylor expansion and defining the Jacobian matrices at time step m.The resulting equation may be expressed as where A, B and J are the Jacobian matrices, and is the variable delta form.In order to reduce the computational effort in the numerical solution of Equation ( 14), the approximate factorization technique (Beam and Warming, 1978) is introduced.Besides, to reach second order spatial accuracy, the scheme uses a central discretization for both convection and diffusion terms; however, to circumvent the odd-even decoupling problem (Patankar, 1980), artificial dissipation terms are added externally.The final form of the equation is: where I being the identity matrix, i  and e  the artifi- cial dissipation coefficients and x  and y  the nodal distance along the corresponding coordinate axes.Following an ADI (alternating direction implicitly) concept, Equation ( 15) is solved as a sequence of two one-dimensional problems Finally, it is important to mention that the artificial dissipation coefficients i  and e  are calibrated numerically and the values found to be adequate were 0.001 i   and 0.003 e   .The complete procedure to control the magnitude of the artificial smoothing terms is described by Zdanski et al. (2004).

Validation/Verification of the Numerical Solution
This section comprises the validation of the numerical method applied.The flow of a binary mixture in a two-dimensional plane channel was the benchmark problem selected, i.e., the turbulent concentration profiles and the Sherwood number distribution were confronted with literature data (for three meshes).Firstly, in Fig. 1 the dimensions of the geometry are presented, i.e., the channel height h=0.011 m and the channel length L=60 h.The computational meshes adopted in the simulations are uniform with 201 × 31, 301 × 35 and 401 × 41 grid points in the streamwise and cross-sectional directions, respectively.Boundary conditions at the inlet plane are: uniform velocity, temperature and concentration profiles, with the Reynolds number based on the hydraulic diameter equal to 27.480.For this test, the following values were applied: u i = 18.68 m/s, C i = 0.00788 kg/m 3 and i  = 0.02 (u i ) 2 /2, with u i , i C and i  being the inlet velocity, concen- tration and turbulent kinetic energy, respectively.At the exit section, parabolic conditions (null derivatives) are specified for all the variables, except for pressure, whose variation is considered to be linear.At the solid walls, a constant mass flux of the chemical species (water vapor) is imposed, " w No = 0.025 kg/m 2 s.   wood number (Sh) expresses the local convective mass transfer coefficient (h m ), i.e., where C m is the bulk concentration of the chemical species, C w is the concentration at the channel walls, h m is the local convective mass transfer coefficient, D h is the hydraulic diameter and u m is the average velocity in the channel cross section.
Re Dh and Sc being the Reynolds and molecular Schmidt numbers.Finally, it is important to mention that the numerical results obtained with the present scheme (see Figures 2 and 3) are clearly mesh independent.

Flow in a Wood Dry Kiln
A typical convective wood dry kiln comprises various channels that are formed between the wood stacks.The previous numerical results of Possamai (2013) and Zdanski et al. (2015), assessing the aerodynamic effects in a dry kiln, have shown that flow inside the channels can be highly non-uniform.Besides, according to the same authors, flow stagnation regions and recirculation zones are typically found in a dry kiln device.Therefore, this section aims at studying the effects of the aforementioned flow characteristics on the drying rate of the wood stack (mass flux in the channel walls).
The geometry adopted in the simulations, along with the flow configuration at the channel entrance section, is shown in Figure 4.A parametric study was performed, i.e., the inlet flow velocity and the incidence angle (Â) were varied to assess the influences on mass transfer at the channel walls.In the present test the boundary conditions were specified similarly to the validation case, except the concentration at the walls that was set to a constant value, C w = 0.025 kg/m 3 , and the inlet flow velocity was varied.The computational mesh adopted in the simulations is uniform with 201 x 31 grid points, the channel dimensions being equal to the validation case.Firstly, in order to study the effects of distinct inlet flow velocities on mass flux at the channel walls, parallel flow was considered (Â=0).The flow velocities tested were u i = 6, 8, 10 and 12 m/s, values typically found in a wood dry kiln (Zdanski et al., 2015).The main results are presented in Figure 5, with the convective mass transfer coefficient represented along the channel walls (see the definition of Equation ( 23)).The analysis of the results clearly indicates that mass flux at the walls is greatly enhanced as the inlet flow velocity varies from 6 to 12 (average increase around 95%).The mass transfer of water vapor is higher at the entrance flow region, reducing to a plateau in the fully developed region (this behaviour is similar to a convective heat transfer phenomenon).Therefore, the main conclusion of the results is that distinct inlet flow velocities affects drastically the mass transfer at the channel walls (and probably the drying rate of the wood stacks).The results of Figure 6 show the convective coefficient for non-vanishing incidence angle at the channel entrance.It is important to mention that the aforementioned flow topology typically occurs inside a wood dry kiln (Possamai, 2013;Zdanski et al., 2015), being worthy of investigation.The entrance velocity u i was set variable with respect to the incidence angle in order to produce an uniform mass flow of the mixture in the channel for the cases tested (Â=0°,10°, 20°, 30° and 40°), i.e., A m udA cte The results of Figure 6 clearly indicate that mass transfer of the chemical species (water vapor) at the channel walls was enhanced in the entrance flow region (x/h<25) as the flow incidence angle increased.The physical scenario is dominated by the stagnation region formed at the channel upper wall, that leads to enhancement in the convective mass transfer process.Otherwise, in the fully developed flow region (x/h>25) the convective coefficient is basically constant (an expected result due to the occurrence of similar velocity profiles in this region).
Finally, in order to conclude this section Figure 7 shows the results for the average Sherwood number along the channel length for Â=10° and 40° and for all velocities tested, i.e., The results show basically a linear variation between Sherwood and Reynolds numbers.The convective mass transfer enhancement was also observed when the flow presented an incidence angle, this behaviour being related to the stagnation region formed at the upper wall.

Flow in Channels with Sudden Expansion
The results of Zdanski et al (2015) show that recirculation regions are formed inside the channels of wood stacks in a dry kiln.Aiming to study the effects of these regions on mass transfer at the walls, the flow in a channel with sudden expansion was idealized.For this test, the boundary conditions adopted are equal to the validation case, except the inlet flow velocity was varied.The geometry is presented in Figure 8, with the channel height h=0.011 m and the channel length L=30 h.The step height s was set variable in a parametric study and the computational Brazilian Journal of Chemical Engineering Vol. 33, No. 02, pp. 287 -296, April -June, 2016 mesh uniform with 901 × 31 grid points in the x and y directions, respectively.The velocity at the inlet section was set as follows: u i = 16 m/s for the channel with expansion ratio ER = (h-s)/(h) = 1/2; u i = 12 m/s for ER = 2/3 and u i = 10 m/s for ER = 5/6.Finally, it is important to mention that u i was set variable with respect to the expansion ratio in order to produce an uniform mass flow of the mixture in a channel for the cases tested (see Equation ( 26)).The purpose of this section was to study the effects of the expansion ratio on the mass flux of the chemical species (water vapor) at the solid walls.Figure 9 shows the results for the convective mass transfer coefficient at the lower and upper channel walls for an expansion ratio ER = 1/2.The convective coefficient at the upper wall presents behaviour similar to the plane channel (the maximum value occurs at the entry region).Otherwise, the recirculation zone formed at the stepped wall produces a distinct behaviour on the convective coefficient, i.e., a maximum value near the reattachment point (similar behaviour was obtained for heat transfer results in a back step problem with the effects of turbulence promoters -see Zdanski et al., 2015) followed by a sharp decrease towards the exit section.The effects of the expansion ratio (ER) on mass flux at the channel walls are presented in Figures 11  and 12.The results clearly indicate that the mass fluxes at both walls are reduced in the entry flow region (x/h < 15) with increasing expansion ratio ER.This behaviour may be explained due to the occurrence of lower velocities in this region since the flow passage area was enhanced.Otherwise, the point of maximum value at the lower wall (see Figure 12) displaces towards the entrance section due to reduction in the recirculation region as the expansion ratio was increased (for the cases simulated x/s = 4, but as the step height s decreases, the vortex length also reduces -see Figures 13 and 14).Finally, in the flow region towards the exit section (x/h > 15) the results for the convective coefficient at the lower and upper walls are similar (with minor differences) for all expansion ratios (this is an expected result since the average velocities inside the channels are equal in this region).

CONCLUSIONS
In this paper a numerical method is presented for solution of the turbulent flow of a binary mixture inside channels.The numerical results obtained indicate that the scheme reproduced satisfactorily the literature data for both concentration profiles and Sherwood number.Furthermore, the numerical solution obtained was clearly mesh independent.
Besides, in the physical analysis performed, the main conclusions are: (i) the parametric study reveals that distinct inlet flow velocities produce drastic effects on the convective mass transfer coefficient.The analysis of the results clearly indicates that mass flux at the walls is highly enhanced as the inlet flow velocity varies from 6 to 12 (an average increase around 95%).Furthermore, the wall mass flux of the chemical species (water vapor) was enhanced in the entrance flow region (x/h<25) as the flow incidence angle increases.The physical scenario is dominated by the stagnation region formed at the channel upper wall, that leads to enhancement of the convective mass transfer process.(ii) The flow in a sudden expansion reveals that the convective coefficient is very sensitive to the expansion ratio ER.In the entry flow region the mass flux reduces as the expansion ratio increases.Furthermore, the maximum value for mass flux near the reattachment point displaces towards the vertical wall as a consequence of the reduction of the vortex zone.
model.In the present model, two differential equations are solved for the turbulence kinetic energy ( ) and its dissipation rate (  ), i.e.,

Figure 1 :
Figure 1: Illustration of the two-dimensional channel with the main dimensions.

Figure 2
Figure 2 shows the computed non-dimensional concentration profiles at the channel exit section (L/h = 60), corresponding to the fully developed flow region.As can be seen in the figure, the results clearly indicate that the numerical scheme reproduces satisfactorily the literature concentration profile for the three meshes tested.The benchmark solution in this case is the classical correlation according to Kader (1981), that was developed to fit the experimental data for fully developed flow regions in plane channels and circular ducts.The non-dimensional concentration is defined as C + =(C-C w )/C, * where C * = " w No /u * and u * = 0.5 0.25 C   (u * being the shear velocity), and the parameter y + (non-dimensional distance from the wall) is given by * y u y     .

Figure 2 :
Figure 2: Comparisons between analytical correlation and computations for the concentration profile in the fully developed flow region.Furthermore, Figure3presents the Sherwood number distribution along the channel walls.The Sher-

Figure 3 :
Figure 3: Comparisons between the empirical correlation and computations for the Sherwood number.The results presented in Figure3indicate a good agreement between the numerical solution and the empirical correlation in the fully developed flow region (x/h  25), the error being around 6%.It is important to mention that the Sherwood number for the fully developed flow region is obtained from the empirical correlation of Dittus-Boelter(Incropera et  al., 2007), i.e.,

Figure 4 :
Figure 4: Illustration of the two-dimensional channel with variable flow incidence angle.
Convective mass transfer coefficient for distinct incidence angles (Â).

Figure 8 :
Figure 8: Illustration of the two-dimensional channel with sudden expansion section.

Figure 9 :
Figure 9: Convective mass transfer coefficient for the channel with expansion ratio ER = 1/2.Aiming to clarify the physics of the problem, Figure 10 presents the streamlines (flow topology) with the vortex region clearly indentified.It is important

Figure 11 :Figure 12 :
Figure 11: Convective mass transfer coefficient for the upper wall in the channel with variable expansion ratios (ER).
transfer coefficient along the walls [m/s] I identity matrix J Jacobian matrices