Unsteady RANS and detached eddy simulation of the multiphase flow in a co-current spray drying☆

Jolius Gimbun ,Noor Intan Shafinas Muhammad ,Woon Phui Law

1 Centre of Excellence for Advanced Research in Fluid Flow,Universiti Malaysia Pahang,Pahang 26300,Malaysia

2 Faculty of Chemical and Natural Resources Engineering,Universiti Malaysia Pahang,Pahang 26300,Malaysia

3 Faculty of Technology,Universiti Malaysia Pahang,Pahang 26300,Malaysia

Keywords:Drying Turbulence Two-phase flow CFD Detached eddy simulation Modelling strategy

A B S T R A C T A detached eddy simulation(DES)and a k-ε-based Reynolds-averaged Navier–Stokes(RANS)calculation on the co-current spray drying chamber is presented.The DES used here is based on the Spalart–Allmaras(SA)turbulence model,whereas the standard k-ε(SKE)was considered here for comparison purposes.Predictions of the mean axial velocity,temperature and humidity profile have been evaluated and compared with experimental measurements.The effects of the turbulence model on the predictions of the mean axial velocity,temperature and the humidity profile are most noticeable in the(highly anisotropic)spraying region.The findings suggest that DES provide a more accurate prediction(with error less than 5%)of the flow field in a spray drying chamber compared with RANS-based k-ε models.The DES simulation also confirmed the presence of anisotropic turbulent flow in the spray dryer from the analysis of the velocity component fluctuations and turbulent structure as illustrated by the Q-criterion.

?2015 The Chemical Industry and Engineering Society of China,and Chemical Industry Press.All rights reserved.

1.Introduction

Spray drying is a dehydration process to convert liquid feed materials into dry powder forms through a hot gas medium.Spray drying is widely used to produce foods,pharmaceutical products and other products such as fertilizers,detergent soap and dyestuffs.

The detailed hydrodynamics of the spray dryer chamber has been studied extensively both experimentally and numerically by several researchers such as Kieviet[1];Kieviet and Kerkhorf[2];Anandharamakrishnan et al.[3];Southwell and Langrish[4];Langrish and Zbincinski[5];Zbicinski et al.[6];Harvie et al.[7];Huang et al.[8].Most of the previous work reported extensive comparison between experimental measurement and computational fluid dynamics(CFD)prediction.Modelling of gas–solid flow in a co-current spray dryer is challenging due to the presence of turbulence,two-phase interactions,and heat and mass transfer.Simulations are often performed using a combination of a simpler two-way coupling gas-solid model and Reynolds-averaged Navier–Stokes(RANS)based turbulence model.Among the RANS-based turbulence models available in commercial FLUENT code,the standard k-ε(SKE)model is the most popular due to its robustness,lower computational demand and ability to give a reasonably accurate prediction.SKE performs well for simple flows and seems to give a fair prediction of the multiphase flow inside the drying chamber.However,there is still a discrepancy on the prediction of gas temperature,axial velocity and the humidity profile,especially in the highly anisotropic spraying region.An advantage of detached eddy simulation(DES)in predicting the flow field in the spray region was successfully demonstrated in this work,whereby the contour plot from DES simulation differs markedly with those from SKE.An accurate prediction of temperature,velocity and humidity profile inside the drying chamber is important,as this region plays an important role in the drying process.It is,therefore,interesting to investigate the capability of various modelling approaches to predict the flow field inside a drying chamber.

The multiphase turbulent flow inside the drying chamber requires a better turbulence model such as the DES.Therefore,this work aims to evaluate the performance of DES in predicting the flow field inside a co-current spray dryer.DES model is a relatively new development in turbulence modelling belongs to a hybrid turbulence model,which blends large eddy simulation(LES)away from the boundary layer and RANS near the wall.This model was introduced by Spalart et al.[9]in an effort to reduce the overall computational effort of LES modelling by allowing a coarser grid within the boundary layers.The DES employed for the turbulence modelling in this work is based on the Spalart–Allmaras(SA)model and has never been previously used for modelling of spray drying.Unlike the SKE,the DES does not suffer from the assumption of isotropic eddy viscosity.Since turbulence flow is anisotropic in nature,thus DES should provide a better prediction of turbulent flow in drying chamber.

2.Spray Dryer Geometry

The spray dryer geometry is shown in Fig.1,with the pressure nozzle atomiser located 229 mm below the top of the chamber,and the drying air enters through an annulus similar to the one studied by Kieviet[1].The air outlet pipe is mounted at the cone centre and is connected to the cyclone to separate the particles from the gas stream.In this work,GAMBIT was used to prepare a three-dimensional computational grid of a co-current spray dryer as illustrated in Fig.1.Predictions from CFD simulations were compared with the laser Doppler anemometer(LDA)measurement by Kieviet[1]at various positions in a spray drying chamber.Data from the CFD simulation were taken as a statistical average up to 1000 time steps(10 s of real time)after a pseudo-steady condition was achieved.Details of the CFD setup are outlined in Table 1.The pressure atomiser model in FLUENT was adopted with spray angle of 76°.The Rosin–Rammler distribution was used to model the particle size distribution using 200 particle classes to represent the spray in the range 10 to 138 μm.The Rosin–Rammler model is given by

where Ydis the retained mass fraction of particle,d is the particle diameter,is the mean particle diameter and n is the size distribution parameter.The feed liquid properties were based on an aqueous maltodextrin solution containing 42.5%solids.The feed liquid has a viscosity of 41.9 mPa·s while the dried particles are often made of a hollow sphere with diameter ranging from 10 to 138 μm.The dried particle size depends on the type of the nozzle used,for instance a twin-fluid nozzle produced a small particle about 20 μm,whereas the pressure atomiser produced larger particles with the mean diameter about 80 μm.The CFD approach used in this work is similar to those described in our earlier work[3],except for the grid and turbulence model employed.The dried particles are collected at the bottom of the cone or through the exit pipe.In addition the particles that come in contact with dryer wall are assumed to be trapped,because most wet droplet may stick to the wall on first contact.

Fig.1.Geometry and surface mesh of co-current spray dryer.

Table 1 Operating condition

3.CFD Approach

3.1.Turbulence model

A two-way coupling method was employed in this work,in which the momentum exchange between both the continuous and discrete phases is taken into account.The liquid droplet feed from the atomiser is assumed to behave as discrete spherical particles,in the similar manner to that of solid particles.This assumption is reasonable for spray dryer where the instantaneous drying of evaporative species from droplet is taking place.Moreover,the droplet is very small in size(10 to 138 μm)and hence the issue of droplet deformation which can affect the particle drag coefficient is not an issue.The pressure atomiser used in this work produced a known droplet size ranging from 10 to 138 μm with the size distribution similar to that of Eq.(1).The droplet was assumed as a mixture containing 57.5%of evaporative species(water)and the remaining content is a non-evaporative species(maltodextrin).

The selection of a turbulence model for spray drying simulation is very important.Extensive model like LES is of course an excellent model,but it is still too computationally expensive to run on a personal computer.Relatively new turbulence models such as DES need to be validated further before they can be applied routinely to spray drying simulation.Therefore,the predictive capabilities of SKE and DES on multiphase flow in a spray dryer have been extensively compared in this study.

The k-ε model is a semi-empirical model based on two transport equations i.e.,the turbulent kinetic energy(k)and its dissipation rate(ε).The Kolmogorov–Prandtl expression for the turbulent viscosity which assumes isotropic turbulence intensity is used.The k-ε model constants according to Launder and Spalding[10]were employed.Turbulence is not resolved for the discrete phase,but rather modelled as stochastic effects of particle interactions with eddy[11].Turbulent particle dispersion is considered in the discrete phase model(DPM)as a discrete eddy concept similar to the one used by Anandharamakrishnan et al.[3].The turbulent air flow pattern is assumed to be made up of a collection of randomly directed eddies,each with its own lifetime and size.

The DES employed in this work is based on the SA model[12].The SA one-equation model solves a single partial differential equation for a variable ? which is called the modified turbulent viscosity:

The variable ? is related to the eddy viscosity by

with additional viscous damping function fv1to ensure the eddy viscosity is predicted well in both the log layer and the viscous-affected region.The model includes a destruction term that reduces the turbulent viscosity in the log layer and laminar sub-layer.

The production term,Gv,is modelled as:

S is a scalar measure of the deformation rate tensor which is based on the vorticity magnitude in the SA model.The destruction term is modelled as:

The closure coefficients for the SA model[12]are Cb1=0.1355,Cb2=0.622,σ?=2/3,Cv1=7.1,Cw1=(Cb1/k2)+(1+Cb2)/σ?,Cw2=0.3,Cw3=2.0,k=0.4187.

The destruction term in Eq.(6)in the SA model is proportional toThe eddy viscosity becomes proportional towhen the destruction term is balanced with the production term.The Smagorinsky LES model varies its sub-grid-scale(SGS)turbulent viscosity with the local strain rate,and the grid spacing is described by vSGSα,where Δ=max(Δx,Δy,Δz).The SA model will act like a LES model if d is replaced with Δ in the destruction term.To exhibit both RANS and LES behaviour,d in the SA model is replaced by:

where Cdesis a constant with a value of 0.65.Then the distance to the closest wall d in the SA model is replaced with the new length scale~d to obtain the DES.The purpose of using this new length is that in boundary layers where Δ by far exceeds d,the standard SA model applies since=d.Away from walls where=CdesΔ,the model turns into a simple one equation SGS model,close to Smagorinsky's in the sense that both make the mixing length proportional to Δ.The Smagorinsky model is the standard eddy viscosity model for LES.On the other hand,this approach retains the full sensitivity of RANS model predictions in the boundary layer.This model has not yet been applied to predict spray drying flows.Applying DES and assessing its performance in relation to experimental data and other turbulence modelling approaches is the main objective of the current study.

3.2.Discrete phase model

The particles-fluid interaction were modelled using a DPM,and two-way coupling was considered in order to enable the prediction of simultaneous heat and mass transfer during the drying process[3,13].The combined Eulerian and Lagrangian model were used to obtain the particle trajectories by solving the force balance equation as follows:

where ugis the fluid phase velocity,upis the particle velocity,ρpis the particle density,ρgis the gas density and g is the gravity.The particle Reynolds number,Re,is given by:

where μ is the fluid viscosity.The drag coefficient,CD,was calculated according to the Morsi–Alexander empirical drag model[14]as follows:

According to Bagchi and Balachandar[15],turbulence does not have a systematic and substantial effect on the mean drag.Therefore,the effect of turbulence on drag is not considered throughout this work.Moreover,the particle is very dilute(less than 1%)in the case of spray drying to affect the continuous phase flow.

3.3.Modelling of heat and mass transfer

The heat and mass transfer between the particles and the hot gas was calculated in the similar manner to Li et al.[16]as follows:

where hfgis the latent heat of vaporization and dmp/dt is the rate of evaporation.The mass transfer between the gas phase and droplet is given by

where kcis the mass transfer coefficient obtained from the Nusselt and Sherwood correlation which is solved by the CFD code.The droplet boiling model is applied to predict the convective boiling of a discrete phase droplet when the temperature of the droplet reached the boiling point while the evaporative species still exists.The boiling rate equation is given by

where cp,gand kgare the gas heat capacity and thermal conductivity,respectively.

3.4.Grid dependent analysis

Study of the grid dependence in CFD calculations of flow field inside the spray dryer was performed in order to find out the minimum mesh density that could yield the acceptable estimations with respect to the experimental measurements.Three different grids(coarse:185 k cells;intermediate:420 k cells;fine:786 k cells)were used to examine the suitability of mesh in this work.These grids consist of hexahedral and tetrahedral meshes.DES turbulent model with unsteady solver was employed for the grid assessment.CFD simulation in this work was performed using six units of HP Z220 workstation with a quad core processor(Xeon 3.2 GHz E3-1225)and 8 Gigabytes of RAM.The CPU time for the coarse grid is below 0.4 s per iteration,whereas the intermediate and fine grids need 0.9 and 1.6 s per iteration,respectively.The results from these three grids were compared with the experimental data from Kieviet[1].Fig.2 shows the axial velocity profile obtained from different grid density.Generally,predictions from these three mesh densities are in good agreement with the experimental data.However,simulation by using coarse mesh(185 k)failed to resolve the double peak flow feature at vertical position Z=1.0 m from the nozzle due to the circular injection of heated air.Both the intermediate(420 k)and fine meshes(786 k)resolved the double peak features accurately.Significant enhancement on the accuracy of prediction was observed when the mesh density was increased from 185 k to higher mesh densities.Hence,the 420 k grid was selected for the remaining of this work to minimize the computational effort.

Fig.2.Comparison of axial velocity between different grid densities with experimental measurement by Kieviet[1].

3.5.Steady and unsteady solver

Most of the CFD studies on spray drying process[1,3,7,8,17]were performed by using steady solver.However,the unsteady solver represents the real measurement better and should be able to produce a more accurate result[18].Similar findings are also reported by Lian and Merkle[18],who found the time-averaged unsteady simulation produced a more accurate prediction on the wall heat flux in a combustion chamber that by steady simulation procedure.Experimental measurement is often taken as time averaged quantities and the unsteady solver mimics this situation much better.The SKE turbulent model was employed to simulate the spray dryer as did in studies by Kieviet[1].The velocity at the point X,Y,Z=0,0,0.3 m was monitored and the time averaging was started only once the velocity at this monitoring point is no longer fluctuating.The time averaging was set for up to 1000 time steps to ensure statistical validity of the data.

Fig.3 shows the comparison of axial velocity predictions from both steady and unsteady solver at various radial positions in the spray dryer.Predictions by both solvers showed good agreement with experimental measurements.However,the prediction from the unsteady solver is much closer to experimental measurements compared to the steady solver,especially at the vertical position Z=1.0 m.This is due to the fact that the experimental measurement was performed by time averaging the instantaneous velocity in the similar way,the unsteady solver was performed.Furthermore,turbulent flow inside the spray dryer is better resolved by using URANS(unsteady solver)than the RANS(steady solver)due to the inherent nature of turbulence.Hence,unsteady solver was employed for the remainder of this work.

Fig.3.Comparison of axial velocity between steady and unsteady solver with experimental measurement by Kieviet[1].

4.Results and Discussion

4.1.Temperature profile

Fig.4 shows the temperature profile versus radial position at the various vertical positions(Z=0.2 m and 1.0 m)of the chamber.The predicted temperature profiles using SKE and DES turbulence models were compared with the data from experimental work[1].Fig.4 showed high temperature fluctuation at the centre region of the spray dryer where the hot air is injected.Closer to the nozzle,the cold droplets made contact with the hot gas with a simultaneous mass and heat transfer activity.As a result,higher temperature fluctuation occurs closer to the nozzle(Z=0.2 m)while lower temperature fluctuations downwards in the spray drying chamber,i.e.Z=1.0 m.The predicted temperature profiles by SKE and DES turbulence models are in good agreement with the experimental measurement by Kieviet[1].Among the turbulence models tested,temperature prediction from the DES model provides a better agreement with the experimental measurement.This may be attributed to the better representation of turbulence by DES by employing a RANS model closer to the boundary layer and the LES model in the bulk region.Fig.5 shows the contour plot of mean axial velocity and temperature.Around the centre region the double peak features can be observed due to the circular inlet for hot gas at the top of the dryer used in this work.The peak temperature of the hot gas inlet also coincides with the hot gas velocity peaks.Temperature and gas velocity decreases further down the drying chamber as the energy was absorbed by the droplet in the form of latent heat of vaporization.

Fig.4.Comparison of temperature between different turbulence models with experimental measurement by Kieviet[1].

Fig.5.Predicted velocity magnitude and temperature contour inside the drying chamber.

4.2.Axial velocity profile

Predicted axial velocity profiles at various positions in the drying chamber is shown in Fig.6.Predictions using both DES and SKE models are in good agreement with the experimental measurement[1].At all vertical positions,predictions from both models show minimal differences except the peaks for the DES model are much higher than those of the SKE model.The differences between both turbulence models tested in this work for the prediction of axial velocity is minimal,with all models capable of predicting the velocity profile very well.This is due to the absence of swirling flow in the chamber,hence it is not critical to use a sophisticated turbulence model to predict the velocity profile in a co-current spray dryer.

Fig.6.Comparison of axial velocity between different turbulence models with experimental measurement by Kieviet[1].

4.3.Humidity profile

Fig.7.Comparison of humidity between different turbulence models with experimental measurement by Kieviet[1].

Fig.7 shows the predicted gas humidity profile at a different vertical distance from the nozzle(Z=0.2 and 1.0 m).The lowest gas humidity in the centre region(?0.2 m ＜R ＜0.2 m)of the spray dryer due to the circular spraying condition of the feed material is predicted correctly using the CFD simulation.The predicted humidity profiles for both turbulence models were in good agreement with Kieviet's measurement[1]throughout the drying chamber.Although,prediction of the SKE model is not as good as the DES model,especially away from the nozzle,e.g.at Z=1.0 m.The prediction error from the DES model is around 5%,probably at the same magnitude of the experimental measurement uncertainties using the micro separator[1].The poor prediction from the SKE model may be attributed by the poor prediction of the temperature profile at Z=1.0 m which in turn affects the mass transfer(evaporation)and heat transfer between particles and the hot air.DES model in another hand has excellent predictions on temperature profile and hence better prediction of humidity profile.There are some minor discrepancies from the DES model predictions at Z=1.0 m,however,it is very difficult to predict humidity profiles accurately,because this CFD model considered the evaporation of moisture from the surface drops to be at a constant drying rate.Droplets may not always be in a constant drying rate regime,especially towards the end of drying.Hence,inter particle diffusion and water desorption factors may also be important to predict humidity.Fig.8 shows the predicted mass fraction of water contour inside the drying chamber.The bulk region of the chamber has about 3.8%mass fraction of water in agreement with the measurement by Kieviet[1]who reported the bulk region of the drying chamber has almost constant temperature and humidity,in exception of the centre region.

Fig.8.Predicted instantaneous humidity contour inside the drying chamber.

4.4.Fluctuating velocity components and turbulent flow structure

Fig.9 shows a contour plot of the predicted fluctuating velocity components in the drying chamber.The axial velocity fluctuation is evidently stronger than both the tangential and radial velocity components.It is,therefore,confirming the presence of anisotropic turbulence in the spray dryer.In most cases,validation of spray drying CFD simulation is presented only for a limited position closer to the gas inlet position,whereby the velocity,temperature and humidity for both DES and RANS SKE simulation does not differ appreciably.In fact,the resultant double peak temperature and velocity feature extends much longer towards the exit pipe.This feature has not been captured by the RANS model.Previous work using RANS SKE also fails to predict the true extent of the temperature and velocity profile for the same case, e.g.Mezhericher et al.[19].

The structure of the turbulence in the drying chamber can be visualised by iso-surfaces of the Q criterion,which is a scalar quantity defined as Q=0.5(Ω2–S2),where Ω is the vorticity magnitude and S is the mean strain rate.Fig.10 shows the iso-surface plot for Q=80,which suggest that the turbulent flow is highly three-dimensional.The turbulent flow is stronger in the centre region of the chamber where the droplets spray and hot gas stream is introduced.The position of the iso-surface plot also coincides with the region where higher velocity fluctuation was observed.

Fig.9.Predicted RMS velocity components.

Fig.10.The structure of turbulence in spray dryer illustrated by iso-surfaces of the Q criterion,Q=80,coloured by velocity magnitude.

5.Conclusions

A detached eddy simulation and an unsteady RANS modelling of gassolid flow in a three-dimensional co-current spray dryer have been simulated.The results suggest that a more accurate prediction of mean velocity,temperature and humidity profile can be obtained using unsteady simulation combined with the DES turbulence model.The DES approach gives a more accurate prediction(with error less than 5%)of temperature and humidity profile in a co-current spray drying especially away from the nozzle region.The DES simulation further confirms the presence of anisotropic turbulence in the spray dryer,hence justifying the demands for a better turbulence model such as DES or LES.The CFD model in this work may be used to further optimise the hydrodynamics in the spray dryer and hence improving product quality.

Nomenclature

Apparticle surface area,m2

a1constant for Morsi and Alexander's drag coefficient

a2constant for Morsi and Alexander's drag coefficient

a3constant for Morsi and Alexander's drag coefficient

Cb1constant of production term

Cb2constant of Eq.(3)

Cgmoisture concentration in the bulk gas,mol·m?3

Cpmoisture concentration at the droplet surface,mol·m?3

Cv1constant of viscous damping function

Cw1constant of destruction term

Cw2constant of Eq.(6)

Cw3constant of Eq.(6)

cp,pspecific heat of particle,J·kg?1·K?1

cp,gspecific heat of gas,J·kg?1·K?1

d distance from wall,m

fv1viscous damping function

fv2turbulence damping function

fwturbulence damping function near wall

Gvproduction term of turbulent viscosity

g gravity acceleration,m·s?2

h heat transfer coefficient,W·m2·K?1

hfgspecific latent heat,J·kg?1

k turbulent kinetic energy,kg·m2·s?2

kcmass transfer coefficient,m·s?1

kgthermal conductivity of the gas,W·m?1·K?1

mpmass of the particle,kg

Re Reynolds number

r dimensionless value

S scalar measure of the deformation tensor

Tggas temperature,K

Tpparticle temperature,K

t time,s

ugfluid phase velocity,m·s?1

uivelocity in i direction,m·s?1

ujvelocity in j direction,m·s?1

upparticle velocity,m·s?1

v molecular kinematic viscosity,m2·s?1

xidistance in i direction,m

xjdistance in j direction,m

xkdistance in k direction,m

Yvdestruction term of turbulent viscosity

Z vertical distance from top of the drying chamber,m

ε turbulent dissipation rate,m2·s?3

ρggas density,kg·m?3

ρpparticle density,kg·m?3

σ?constant for characteristic stress

χ constant of viscous damping function