Simulation of a single spherical droplet drying considering shrinking and temperature profiles: first drying stage

Abstract

A model including simultaneous droplet heating and water evaporation is proposed to simulate temperature, shrinking and mass profiles of a spherical droplet subjected to convective drying, being valid for the first drying stage. Experimental data on drying skim milk and colloidal silica obtained in the literature were used for validation, but there is no restriction in the model that prevents it from being suitable for other materials. There were not significant differences observed concerning to the droplet components (dissolved or insoluble materials). The initial heating time of the particle upon reaching the constant temperature is relatively short ($\Delta t\approx7s$) for both simulated materials and water evaporation during the first drying stage occurs mostly at the wet bulb temperature of the air. Discrepancy between simulated and experimental values did not exceed $9\%$ for skim milk and $7\%$ for colloidal silica in this first stage, indicating good applicability of the model. Considering the applicability of the model in a more generic way, Whitaker correlation evaluated at the film temperature showed better results. Finally, the small discrepancy found is discussed and some improvements are proposed.

Key words
Critical moisture; differential-algebraic equations; first drying stage; single droplet drying; spray drying

Introduction

Drying can be defined as an operation that aims to remove the volatile liquid from a solid medium through a process involving simultaneous heat and mass transfer. In general, drying occurs in a convective manner by the effective contact between solid particles (including droplets of suspension) and a low humidity gas which is usually heated. For food materials, moisture must be reduced to less than 10% to avoid microbial growth that can negatively affect the quality of the dried product (Yilmaz et al. 2019, González-Toxqui et al. 2020). In spray drying - commonly used in the drying of food, dairy products, inorganic chemicals, ceramic powders and nanoparticles - atomization is the most important step that defines the properties of liquid droplets after spraying, resulting in several forms of the dry product from fine powders to granules (Campos et al. 2021).

Once there are simultaneously numerous droplets during the spray drying process, it is an extremely difficult task to accurately follow the evolution during the drying of these particles inside the dryer (Abdullahi et al. 2020). In the inner of spray dryer chambers, in continuous processing, there are particles exposed to the drying air with different residence times and, therefore, with different humidity at the same time. The kinetics of the drying process of a population of particulate materials is based on the history of humidity and temperature of individual particles (Keshani et al. 2015). Many scientific studies have investigated single droplet drying (SDD) on an experimental basis. These studies allow analysis of the effect of parameters during drying such as temperature, size, morphology and composition (Schutyser et al. 2018). Considering the drying behavior of a single particle, average humidity and temperature values can be estimated for the powder product concerning particles of different residence times inside the dryer through a population balance (Handscomb et al. 2009, Birchal & Passos 2005).

Since there is difficulty in carrying out single droplet drying experiments, many studies are based on literature experimental data. In one of the most classic research projects on SDD experiments (Nešić & Vodnik 1991), temperature and mass data are available for skim milk and colloidal silica particles, making this experiment data widely used for validation of simulated models. Single droplet drying studies, however, demand more in-depth analysis to understand the properties of particles with a high content of solids regarding to its crust formation during drying, especially when the solution is multicomponent (Schutyser et al. 2018). The shrinkage behavior of skim milk droplets with 50% solids as a function of the drying air temperature was examined by Fu et al. (2013). With high initial solids concentrations, the shrinkage of these droplets is directly influenced by the fraction of solids since the beginning of drying, affecting the structural formation of the particle.

The maximum structural stability of the particle exposed to spray drying is spherical. However, changing the drying process parameters, which are related to the mass and heat transfer of the droplet, can result in different particle morphologies (Vicente et al. 2013). Different morphologies depend on how shrinkage occurred, which is influenced by material properties, drying rate and initial diameter. If drying is rapid, solids tend to accumulate in the outer layer of the droplet, where the particle exhibits less shrinkage as the crust formed is practically solid (Chew et al. 2015). In the classic experimental tests carried out by Nešić & Vodnik (1991), the suspended droplets were practically spherical and the authors estimated the error of treating them as spheres as less than $1\%$. The morphology of magnetic iron oxide particles coated using spray drying was investigated and through image and X-ray analysis, the verified particles were approximately spherical in shape while the surface of these particles is composed only of the coating material, as expected (Donadel et al. 2009). Analyzing single droplet drying experiments, the development of the particle morphology of whey protein-maltodextrin mixtures may be related to its rheological behavior, and the developed models describe that viscosity increases as a function of concentration. In addition, morphology and skin formation are apparently dependent on drying time (Both et al. 2019).

Therefore, the knowledge of how water migrates during the drying of the droplets is extremely important for understanding the formation of the particle, its morphology and composition, which significantly affects the quality of the final product. Through simulation, it is possible to have much more accurate analyzes. Differently from what is generally discussed in the bibliography, using more accurate numerical methods, it is not necessary to make some considerations for simplicity, given that the methodology proposed in this work allows representing real drying behavior in a more reliable way. The aim of this study is to investigate the first stage of the spherical particle drying process by proposing and simulating a dynamic mathematical model including heating, water evaporation and particle shrinkage, simultaneously. Model validation was performed by Chaves et al. (2019) using experimental skim milk drying data available in the literature. Here, the proposed model is validated for another material, colloidal silica, and the correlations used to estimate the heat and mass transfer coefficients are verified.

MATERIALS AND METHODS

A mathematical model is proposed and for its simulation, integral balances in control volumes (finite volume methods) were performed and a differential-algebraic equations system is proposed. Constants and numeric values of parameters necessary to simulate the model were synthesized by Dalmaz et al. (2007) for skim milk and colloidal silica. Experimental data were taken from Nešić & Vodnik (1991) for the same products to validate the phenomenological model. In the experiments discussed, droplets were suspended in a controlled air current and their weight and temperature were measured (in separate experiments). For each evaporation curve, $10-20$ practically identical droplets were tested.

Initial considerations

Typical spray drying process for a particle can be divided into two stages. In the first one, heat is supplied to the droplet and the moisture is removed at a constant rate. Basically, due to the diffusion mechanism, moisture migrates to the surface of the droplet, evaporating at a rate sufficient to maintain the air saturation around it. In this first stage, the drop volume change corresponds to the amount of evaporated water. When the drop reaches its critical moisture content (d.b.), $X_{cr}$, the formation of a porous crust on its external surface begins. This provides some resistance to the diffusion of water and, consequently, the drying rate becomes decreasing. In this second stage, in a more common consideration, the volume of the particle remains practically constant with the decrease in humidity, however, the thickness of the crust increases, as does the temperature of the particle. This occurs until the humidity inside the particle is equal to that of its surface, which is in equilibrium with that of the drying air. Critical humidity, $X_{cr}$, represents the point at which the surface free water has been evaporated and the water inside the particle’s pores begins to be removed by diffusion mechanisms (Mezhericher et al. 2007, Handscomb et al. 2009).

Mezhericher et al. (2007) suggested that the first drying stage may be subdivided into two: (i) initial heating of the droplet, being fast enough that no significant evaporation occurs during this period and that the droplet radius $R_d$ remains constant. However, Mezhericher et al. (2007) also pointed out that the temperature distribution inside the drop must be considered; (ii) evaporation of water on the surface of the drop after its initial heating. If dimensionless number of Biot is less than 0.1 and dimensionless number of Fourier is more than $0.1$, the global capacitance method can be used. Otherwise, the temperature distribution inside the particle must also be examined. The radius $R_d$ varies with time since the shrinkage rate is proportional to the evaporation rate. These authors disregard the heat transfer by radiation since Lin & Chen (2002) estimated that, at most, this form of heat exchange represents $4\%$ of the total heat transfer, which occurs basically by convection. As the droplet develops a crust it would be more logical to call it a particle (Nešić & Vodnik 1991).

Although Mezhericher et al. (2007) neglected the water evaporation on the droplet surface during its initial heating, it is a fact that the evaporation occurs as soon as the particle encounters the drying air. Thus, even in the initial heating phase, water is evaporated and the droplet reduces its diameter. It is then possible to develop a model that will describe the first stage of drying a spherical particle comprising simultaneously its heating, the water evaporation and its shrinkage.

Model presentation

The proposed model presented is valid to describe the first droplet drying stage, that is, until the particle reaches its critical humidity. The distribution of temperature inside the droplet is considered, as well as the variation of its properties over time. As suggested by Mezhericher et al. (2007), radiation transfer of heat is neglected. A convenient way to express the energy balance is as function of enthalpy (Equation 1) (Maliska 2013). The substantive derivative $\left(D/Dx\right)$ involves the specific mass $\rho$ and specific enthalpy $h$ with velocity $u_r$ (considering only radial velocity), which would encompass mass transfer within the particle (Equation 2).

$$\frac{D\left(\rho h\right)}{Dt}=\dot{q}-𝐝𝐢𝐯\vec{q}+\frac{Dp}{Dt}+\mu \varphi$$ (1)

$$\frac{D\left(\rho h\right)}{Dt}\equiv \frac{\partial \left(\rho h\right)}{\partial t}+\frac{\partial \left(\rho u_{r}h\right)}{\partial r}$$ (2)

There is no internal energy generation ($\dot{q}$ is zero) and the droplet pressure $p$ variation is negligible (so the substantive derivative of pressure is zero). The viscous effects of the fluid are also disregarded ($\mu \varphi$ is also zero). To express heat transfer, $𝐝𝐢𝐯\vec{q}$, Fourier’s law is used in spherical coordinates, considering only the variation of properties with radius $r$. Thus, the differential equation that governs the process is presented in Equation 3 together with the boundary conditions (Equation 4) and initial conditions (Equation 5). Fourier’s law expresses heat flux in terms of droplet thermal conductivity $k_d$ and temperature gradient $\partial T_d/dr$. The boundary conditions consider symmetry at the droplet center, as well as convection and evaporation only at the surface ($r=R_d$).

$$\frac{\partial \left({\rho }_{d}h\right)}{\partial t}=\frac{1}{r^2}\frac{\partial }{\partial r}\left(k_d{r}^2\frac{\partial T_d}{\partial r}\right)-\frac{\partial \left({\rho }_d{u}_{r}h\right)}{\partial r}$$ (3)

$$\{\begin{array}{cc}\hfill h_c(T_{\infty }-T_d)=k_d\frac{\partial T_d}{\partial r}+\lambda \frac{\dot{m_v}}{A_d}\hfill & \hfill forr=R_d\left(t\right)\hfill \\ \hfill \frac{\partial h}{\partial r}=0\hfill & \hfill forr=0\hfill \end{array}$$ (4)

$$\{\begin{array}{cc}\hfill \begin{array}{c}\hfill r=R_{d,0}\hfill \\ \hfill T_d=T_0\hfill \end{array}\hfill & \hfill fort=0\hfill \end{array}$$ (5)

For the development of this model, the radial direction is divided into $N$ equally spaced points, with the first point in the center and the last one on the surface of the particle. Thus, the particle is divided in the radial direction into control volumes, which, except for the outermost, have their average properties calculated at the central point of each volume, $r_i$, as shown in Figure 1a. It can be seen in Figure 1a that the properties of the innermost volume are calculated at the central point of the sphere ($r_1=0$). For the outermost volume, the properties are calculated on the particle surface ($r_N=R_d$), as shown in Figure 1b. The temperature derivatives on the faces of the volumes were equated by finite central differences.

Figure 1
Spherical particle discretization for generic point $i$ (a), and outer surface (b).

The differential Equation 3 is discretized (with $i$ being the discretization index) and integrated into a differential spherical shell volume $dV$ (Equation 6). By considering symmetry around the center, the properties do not vary with the angles $\phi$ and $\theta$, varying only with the radius, so the differential element of integration can be simplified (Equation 7) and the integral is calculated only on the radius with limits of integration from $r_{i-1/2}$ to $r_{i+1/2}$ in order to express the discretized balance for each control volume. The indexes $i-1/2$ and $i+1/2$ indicate, respectively, the internal and external faces of the control volume $i$.

$$\iiint \frac{\partial \left({\rho }_{d}h\right)}{\partial t}dV=\iiint [\frac{1}{r^2}\frac{\partial }{\partial r}\left(k_d{r}^2\frac{\partial T_d}{\partial r}\right)-\frac{\partial \left({\rho }_d{u}_{r}h\right)}{\partial r}]dV$$ (6)

$$dV=r^{2}sin\phi d\phi d\theta dr=4\pi r^{2}dr$$ (7)

The Equation 8 describes the energy balance in each internal volume of the droplet. $H$ represents the enthalpy of the control volume and $\dot{m}$ the mass transfer rate. The sub-index $w$ is for water and $d$ is for droplet,

$$\frac{\mathrm{\partial }{H}_i}{\mathrm{\partial }t}=-4\pi r_{i-\frac{1}{2}}^2{k}_{d,i-\frac{1}{2}}{\frac{\mathrm{\partial }{T}_d}{\mathrm{\partial }r}|}_{i-\frac{1}{2}}+4\pi r_{i+\frac{1}{2}}^2{k}_{d,i+\frac{1}{2}}{\frac{\mathrm{\partial }{T}_d}{\mathrm{\partial }r}|}_{i+\frac{1}{2}}+{\dot{m}}_{w,i-\frac{1}{2}}{h}_{w,i-\frac{1}{2}}-{\dot{m}}_{w,i+\frac{1}{2}}{h}_{w,i+\frac{1}{2}}$$ (8)

By discretizing the particle, it is necessary to include the effect of mass transfer of water within the particle. The last two terms of Equation 8 correspond to energy transport between adjacent discretized volumes due to liquid water transfer. One consideration made is that the interface temperature between adjacent control volumes is calculated as the average of the respective temperatures of each volume (Equations 9 and 10). The radii are calculated in an analogous way (Equations 11 and 12).

$$T_{i-\frac{1}{2}}=\frac{T_i+T_{i-1}}{2}$$ (9)

$$T_{i+\frac{1}{2}}=\frac{T_i+T_{i+1}}{2}$$ (10)

$$r_{i-\frac{1}{2}}=\frac{r_i+r_{i-1}}{2}$$ (11)

$$r_{i+\frac{1}{2}}=\frac{r_i+r_{i+1}}{2}$$ (12)

The use of the enthalpy of the particle instead of its temperature facilitates the resolution of the model, since the moisture content $X$ and, consequently, the mass of each control volume varies with the drying time. With this, the time derivative has the product of several time functions ($V$ - volume, $\rho$ - specific mass, $c_{p,d}$ - specific heat), in addition to temperature ($T$). Furthermore, if the correlation that provides some property is modified, it will not be necessary to re-equate the entire model.

In addition to the $N$ differential equations of energy balances (Equation 8), the model has $N$ algebraic equations that define the enthalpy of each control volume (Equation 13) and a differential equation for determining the variation of the droplet diameter (Equation 32). In this way, there is one differential-algebraic equations system with $2N+1$ equations. The numerical resolution of the DAE system follows a methodology similar to that implemented by Costa Jr et al. (2003), who have highlighted that this type of resolution facilitates algebraic manipulation and reduces the consumption of computational time without loss of information through differentiation. The enthalpies for each control volume are calculated by Equations 13 and 14.

$$H_i={\rho }_{d,i}\left(X\right)V_i{h}_i$$ (13)

$$h_i={\int }_{T_{ref}}^{T_i}{c}_{p,d}(T,X)dT$$ (14)

For the center of the droplet, the heat rate is supplied from $V_2$ to $V_1$ and then the water rate is transferred from $V_1$ to $V_2$ (Equation 15). For spherical shells from $i=2$ to $i=N-1$, the energy balance (Equation 8) is used. At the surface, $i=N$, the energy balance presented by Equation 16 describes that heat enters the droplet by convection and mass evaporates from the outermost spherical layer of the droplet to the gaseous medium being $\lambda$ the latent heat of water vaporization. As the droplet shrinks during the first drying step, its radius $R_d$ varies with time. The sub-index $v$ corresponds to vapor and since the only volatile substance is water and as shrinkage is considered uniform, the mass rate of water transferred at the interfaces ${\dot{m}}_w$ between each of the various internal control volumes is proportional to the mass rate of water evaporated on the surface $\dot{m_v}$ and to each volume in relation to the droplet total volume.

$$\frac{\mathrm{\partial }{H}_1}{\mathrm{\partial }t}=4\pi r_{1+\frac{1}{2}}^2{k}_{d,1+\frac{1}{2}}{\frac{\mathrm{\partial }T}{\mathrm{\partial }r}|}_{1+\frac{1}{2}}-{\dot{m}}_{w,1+\frac{1}{2}}{h}_{w,1+\frac{1}{2}}$$ (15)

$$\frac{\mathrm{\partial }{H}_N}{\mathrm{\partial }t}=-4\pi r_{N-\frac{1}{2}}^2{k}_{d,N-\frac{1}{2}}{\frac{\mathrm{\partial }T}{\mathrm{\partial }r}|}_{N-\frac{1}{2}}+4\pi {R_d}^2\left(t\right)\left[h_c(T_{\mathrm{\infty }}-T_{d,surf})\right]-\lambda \dot{m_v}+{\dot{m}}_{w,N-\frac{1}{2}}{h}_{w,N-\frac{1}{2}}-{\dot{m}}_v{h}_{w,N}$$ (16)

It should be noted that the last two terms of the Equation 16 correspond to the transport of energy in the form of sensible heat from one element to the other due to the water transfer. In addition to the energy loss in the surface by evaporation, mass of liquid water is transferred between adjacent volumes, and energy is carried with it.

There are always errors intrinsic to the process of discretizing continuous equations. The finite difference method used in this study uses the central difference approximation with an error in the order $O\left(\Delta r^2\right)$, where $\Delta r$ is the spatial distance between two discretization points. The mesh of the discretized points varies in size due to the droplet shrinkage, which is considered uniform. The calculations of $r_i$ (equation 17) and $\Delta r$ (Equation 18) must therefore be readjusted at each iteration. The discretization error can be reduced by increasing the number of finite volumes $N$. The mesh convergence test was performed based on mass and temperature profiles, observing which value of $N$ that satisfies the convergence criterion. This type of mesh convergence analysis can allow obtaining an accurate solution with a mesh that is dense enough and does not require excessive use of computing resources.

$$r_i=(i-1)\Delta r=(i-1)\frac{R_d\left(t\right)}{N-1}$$(17)

$$\Delta r=\frac{R_d\left(t\right)}{N-1}$$(18)

The convective heat transfer coefficient $h_c$ is estimated by the correlation of Ranz & Marshall (1952) (Equation 19) where the air properties are evaluated at $T_{\infty }$ noting that the correlation is valid for Reynolds numbers between $0$ and $200$.

$${Nu}_d=2+0.6{Re}_d^{\frac{1}{2}}{Pr}^{\frac{1}{3}}$$ (19)

An alternative correlation to Equation 19 is a correlation proposed by Whitaker (1972) (Equation 20), tested for $0.71<Pr<380$ and $3.5<Re<7.6\times {10}^4$. All properties are evaluated at the air temperature $T_{\infty }$ except ${\mu }_{air,surf}$ which is the air viscosity evaluated at the surface temperature $T_{surf}$ of the spherical particle. Thus, the boundary layer effects are incorporated in the determination of the convective coefficient.

${Nu}_d=2+(0.4{Re}_d^{\frac{1}{2}}+0.06{Re}_d^{\frac{2}{3}}){Pr}^{0.4}{\left(\frac{{\mu }_{air}}{{\mu }_{air,surf}}\right)}^{\frac{1}{4}}$ (20)

The specific heat of the droplet considering it as an ideal mixture of solids and water can be evaluated by Equation 21 where $c$ is the concentration of solids that depends on moisture $X$. The sub-index s refers to the solid material of the droplet and $0$ represents the initial condition of the parameter.

$$c_{p,d}=c_{p,w}(1-c)+c_{p,s}c$$ (21)

$$c=\frac{1}{1+X}$$ (22)

$$X=\frac{m_d(1+X_0)}{m_{d,0}}-1$$ (23)

The thermal conductivity of the droplet $k_d$ can be determined as a combination of the thermal conductivities of the water and the solid (Equation 24) where the void fraction $\delta$ is defined as the ratio between the volume of water in the particle and the total volume of the particle (Equation 25). The specific mass of the droplet ${\rho }_d$ is calculated as a function of moisture content $X$ (Equation 26).

$$k_d=\delta k_w+(1-\delta )k_s$$ (24)

$$\delta =\frac{V_{d,w}}{V_d}$$ (25)

$${\rho }_d=\frac{(1+X){\rho }_{d,s}{\rho }_{d,w}}{({\rho }_{d,w}+X{\rho }_{d,s})}$$ (26)

The droplet mass transfer rate is calculated on the surface area $A_d$ and it is determined using Equation 27 with mass transfer coefficient $h_m$ estimated by Equations 28 or 29. Ranz Marshall and Whitaker correlations for calculating the heat transfer coefficient can be applied to mass transfer problems replacing ${Nu}_d$ for ${Sh}_d$ and $Pr$ for $Sc$.

$$\dot{m_v}=h_m({\rho }_{v,surf}-{\rho }_{v,\infty })A_d$$ (27)

$${Sh}_d=2+0.6{Re}_d^{\frac{1}{2}}{Sc}^{\frac{1}{3}}$$ (28)

$${Sh}_d=2+(0.4{Re}_d^{\frac{1}{2}}+0.06{Re}_d^{\frac{2}{3}}){Sc}^{0.4}{\left(\frac{{\mu }_{air}}{{\mu }_{air,surf}}\right)}^{\frac{1}{4}}$$ (29)

The values of water vapor densities on the surface of the droplet, ${\rho }_{v,surf}$, and in the surrounding fluid, ${\rho }_{v,\infty }$, are calculated as a function of temperature, according to Equations 30 and 31. Drying air characteristics such as velocity $v_{air}$ and absolute humidity $H_{air}$ are important parameters to calculate heat and mass transfer as well as some constants such as molar mass $M$ and the universal gas constant $R$.

$${\rho }_{v,\infty }=\frac{H_{air}{P}_{atm}{M}_{air}}{R(T_{\infty }+273.15)}$$ (30)

$${\rho }_{v,surf}=\frac{M_w{P}_{sat}}{R(T_{d,surf}+273.15)}$$ (31)

The shrinkage rate is calculated by Equation 32 as a function of the evaporation rate at the droplet surface (Equation 27). By integrating the differential equation of droplet shrinkage (Equation 32), it is possible to determine the droplet mass (Equation 33) as a function of its actual diameter and its initial values of mass and diameter.

$$\frac{d\left(d_d\right)}{dt}=-\frac{2}{{\rho }_{d,w}\pi d_d^2}\dot{m_v}$$ (32)

$$m_d=m_{d,0}-\pi {\rho }_{d,w}\frac{(d_{d,0}^3-d_d^3)}{6}$$ (33)

The correlations used for air properties are presented in Equations 34 to 37. The correlations for ${\rho }_{air}$, ${\mu }_{air}$ and $k_{air}$ were taken from Cheong (1983) and $c_{p,air}$ from Sandler (2020).

$${\rho }_{air}=1.293\left(\frac{273.15}{273.15+T}\right)$$ (34)

$${\mu }_{air}=1.720\times {10}^{-5}+4.568\times {10}^{-8}T$$ (35)

$$k_{air}=1.731(0.014+4.296\times {10}^{-5}T)$$ (36)

$$c_{p,air}=969.542+6.801\times {10}^{-2}(T+273.15)+16.569\times {10}^{-5}{(T+273.15)}^2-67.828\times {10}^{-9}{(T+273.15)}^3$$ (37)

The correlation used for the water diffusion coefficient in the air is described by Equation 38. The correlation for calculating the saturation pressure of liquid water is expressed by Equation 39.

$$D_{w,air}=0.220\times {10}^{-4}{\left(\frac{273.15+T}{273.15}\right)}^{1.75}$$ (38)

$$P_{sat}=133.3exp(18.3036-\frac{3816.44}{T_{surf}+229.02})$$ (39)

The set of Equations 8 to 32 form a system of differential-algebraic equations (DAEs) and constitute the model to be simulated in this work. This system has no analytical solution and must be solved numerically. The system of differential-algebraic equations was integrated over time using Matlab’s ode15i function. Dalmaz et al. (2007) presented synthesized tables with the necessary parameters for simulation the model. Experimental data obtained by Nešić & Vodnik (1991) was also used to simulate and validate the single droplet drying of skim milk and colloidal silica in a stream of hot air with constant properties.

Milk is a complex food, with many different molecular species, consisting basically of water, fat, lactose, proteins and minerals. The skim milk droplet is considered to be an ideal mixture of water (volatile component) and dissolved solids. The data presented in Table I details skim milk and water property values that were used in the simulation of the drying model. The skim milk droplet with initial mass fraction of solids of approximately $0.2$ - which corresponds to $X_0=4.0$ - was dried at a temperature of $50$C.

Similar to skim milk, a single droplet of colloidal silica with initial mass fraction of solids of approximately $0.3$ - which corresponds to $X_0=2.333$ (d.b.) was dried at a temperature of $101$°C. However, colloidal silica particle is considered to be an ideal mixture of water and insoluble and dispersed solids. Chemically, silica is silicon dioxide ($Si{O}_2$), a material usually inert and insoluble in most solvents, including water. Colloidal silica suspensions with concentrations between $30\%$ and $50\%$ of solids are usually available on the market. However, at more than $40\%$ concentrations of $Si{O}_2$, colloidal silica turns into a gel structure (Nešić & Vodnik 1991). Data presented in Table II were used in the simulation and validation of the colloidal silica drying model. The parameters for water ${\rho }_{d,w}$, $k_w$, $c_{p,w}$ and $\lambda$ remain the same as those used for the skim milk drying simulation (see Table I).

Table I
Parameters used for drying of skim milk (Dalmaz et al. 2007).

Parameter	Value	Unit
${\rho }_{d,s}$	$1281$	kg/m³
$k_s$	$0.372$	W/(m K)
$c_{p,s}$	$1791$	J/(kg K)
${\rho }_{d,w}$	$1000$	kg/m³
$k_w$	$0.620$	W/(m K)
$c_{p,w}$	$4178$	J/(kg K)
$R_{d,0}$	$8.684\times {10}^{-4}$	m
$X_0$	$4.0$	kg(water)/kg(dry solids)
$X_{cr}$	$1.565$	kg(water)/kg(dry solids)
$\lambda$	$2200\times {10}^3$	J/kg
$v_{air}$	$0.75$	m/s
$H_{air}$	$0.01$	kg/kg
$T_0$	$17$	°C
$T_{\infty }$	$50$	°C

Table II
Parameters used for drying of colloidal silica (Dalmaz et al. 2007).

Parameter	Value	Unit
${\rho }_{d,s}$	$2220$	kg/m³
$k_s$	$1.445$	W/(m K)
$c_{p,s}$	$825$	J/(kg K)
$R_{d,0}$	$9.445\times {10}^{-4}$	m
$X_0$	$2.333$	kg(water)/kg(dry solids)
$X_{cr}$	$1.287$	kg(water)/kg(dry solids)
$v_{air}$	$1.73$	m/s
$H_{air}$	$0.01$	kg/kg
$T_0$	$29.3$	°C
$T_{\infty }$	$101$	°C

Simulation is automatically interrupted when the droplet moisture content reaches the critical value $X_{cr}$ since this work models and simulates only the first single droplet drying stage. To compare simulated results to experimental data, the average ($avg$) temperature of the droplets is calculated in spherical coordinates (Equation 40) based on the temperature distribution within the droplet.

$$T_{d,avg}=\frac{3}{{R_d}^3}{\int }_0^{R_d}{r}^{2}T\left(r\right)dr$$(40)

RESULTS AND DISCUSSIONS

The proposed model presented and validated involves enthalpy balances and the numerical resolution comprises an implicit method of solving differential equations. The use of the enthalpy of the droplet instead of its temperature facilitates the resolution of the model, since several properties other than temperature vary with drying time and any necessary change in correlations and/or parameters does not entail a re-equation of the entire model.

Single droplet drying and shrinkage profiles are presented in this section for skim milk and colloidal silica. Spatial distributions of temperature over time are also discussed. The model was previously validated by Chaves et al. (2019) with experimental data obtained by Nešić & Vodnik (1991) for skim milk single droplet drying using the Ranz & Marshall (1952) correlation to estimate the heat and mass transfer coefficients. In this work, another correlation is analyzed (Whitaker 1972), in addition to validating the proposed model for drying another product, colloidal silica, which differs from milk because it is considered an inorganic product containing insoluble solids. It is proposed, therefore, to observe the drying behavior during the first stage and compare drying of the two materials as well as possible simulation differences obtained by each correlation (Ranz-Marshall and Whitaker).

In a preliminary analysis, it was observed that the two correlations (evaluated at $T_{\infty }$) used were not suitable for the silica drying simulations. Colloidal silica was dried in a relative high temperature drying air, $101$°C. There is a considerable large temperature gradient that the correlations do not consider. Therefore, it was decided to modify the correlations and also evaluate them at the average temperature $\bar{T}$ (also called film temperature) between the drying air temperature $T_{\infty }$ and the droplet surface temperature $T_{surf}$.

The mesh convergence for the numerical solution provided was analyzed and the discretization errors tend to vanish as the mesh size is decreased. The mesh convergence analysis was performed both for mass and temperature profiles, observing which value of $N$ that satisfies the convergence criterion. The relative error in relation to a large number of finite volumes ($N=50$) was used as a convergence criterion, since graphically and numerically the convergence can already be obtained from $N=4$. The errors presented in Figure 2 are of order much smaller than the variables of interest, and it is also important to emphasize that these are all smaller than the uncertainties of the measuring instruments used by Nešić & Vodnik (1991). After analyzing the results, mesh convergence is considered from $N=10$ and in this work $N=15$ was used to have adequate and accurate results and still save simulation time.

Figure 2
Mesh convergence analysis using error criteria for the two variables of interest: average temperature (a) and droplet mass (b).

To investigate the performance of the simulated models, RMSE (root mean squared error) and the coefficient of determination $R^2$ were evaluated. Such indicators are shown in Table III for the temperature variable. The RMSE has one unit (dimension) equal to the dimension of the observed and predicted values. $R^2$, in turn, represents, in percentage terms, the variation in the response that is explained by the model. For all simulations presented in Table III, the Reynolds number values varied between $54$ and $178$, which characterizes laminar flow regimes and provides the applicability of both correlations for all cases studied.

Table III
Performance evaluation (RMSE [C] and ${𝐑}^2$) for droplet average temperature simulations.

Product	Correlation	Evaluated at ${𝐓}_{\infty }$		Evaluated at $\overline{𝐓}$
Skim milk	Ranz-Marshall	$RMSE=1.0787$	$R^2=0.7960$	$RMSE=1.1026$	$R^2=0.7869$
	Whitaker	$RMSE=1.0418$	$R^2=0.8098$	$RMSE=1.0584$	$R^2=0.8036$
Colloidal silica	Ranz-Marshall	$RMSE=1.6834$	$R^2=0.2074$	$RMSE=1.0077$	$R^2=0.7160$
	Whitaker	$RMSE=1.6329$	$R^2=0.2543$	$RMSE=0.9635$	$R^2=0.7403$

The dynamic profiles of the average drying temperature of the skim milk droplet (Figure 3a) and colloidal silica (Figure 3b) presented were for the best conditions presented in Table III, that is, the ones with the smallest errors, namely, Whitaker-$T_{\infty }$ for skim milk and Whitaker-$\bar{T}$ for colloidal silica. Experimental data taken from Nešić & Vodnik (1991) were used to validate the proposed model. Assessing the results in Table III, despite the Whitaker correlation has presented better performances in general, for milk, all simulations were satisfactory and obtained results with high concordance between simulated and experimental data, with $R^2$ values considered high. Whitaker (1972) discusses the choice of the reference temperature. Film temperature as the reference temperature would eliminate the need to include some parameter relating to fluid viscosity in the correlation. However, the author chose to present a more generic correlation that can be applied to different types of fluids, including those more viscous than air. In this present work, for silica, what effectively resulted in an improvement in the model was the use of the film temperature $\bar{T}$, regardless of the correlation used. This is due to the fact that the drying air was at $101$°C while the droplet stayed most of the time at an approximate temperature of $32$°C. This temperature gradient is high and not explained by correlations. The use of $\bar{T}$ indirectly included the influence of the boundary layer.

Figure 3
Time profile of the average temperature of a single droplet drying of skim milk with drying air at 50C (a), of colloidal silica with drying air at 101C (b).

Temperature profiles (Figure 3) show an initial period of rapid heating and then a plateau in which the temperature remains practically constant. This plateau is caused by the water evaporation from the droplet surface as the water molecules absorb energy from the drying air to give them the necessary movement to escape the surface and become a gas. Water evaporation during this stage is called quasi-equilibrium evaporation and takes place at a slightly higher value than the wet bulb temperature, which is the lowest temperature that can be reached by evaporating water. It is observed, in general, that the simulated temperature values describe the experimental data satisfactorily (deviation $<9\%$ for skim milk and $<5\%$ for colloidal silica). It is important to point out that in the experiments by Nešić & Vodnik (1991), the temperature was followed by a micro thermocouple that firstly, when calibrated with evaporation of pure water droplets, had an error between $1$ and $7$C, and after compensation, this error dropped to approximately $0.5$C. The simulation ends in $t\approx 150s$ for skim milk and in $t\approx 39s$ for colloidal silica, when the critical moisture content value is reached, which indicates the end of the first drying stage. For $t>20s$ in milk drying and $t>10s$ in silica drying, the simulated data show that the droplet temperature tends to a constant temperature which corresponds to the wet bulb temperature of the air. Numerical values of this temperature are $23.0$C for drying skim milk at $50$C and $32.6$C for drying colloidal silica at $101$C, both with air absolute humidity of $10$g/kg.

The prediction of the average droplet temperature during the rapid droplet heating is adequate, but the model underestimates the rate of temperature increase when the crust formation begins. This may be due to the inaccuracy in the used value of $X_{cr}$. According to Dalmaz et al. (2007), this value was obtained from the experimental curve of drying the droplet in a gas stream. It can be seen in Figure 3, an increase in the particle temperature at $t>80s$ for milk and at $t>10s$ for silica, indicating the beginning of the second drying stage at a decreasing rate. Birchal & Passos (2005) specified the value of $X_{cr}$ in the spray drying of whole milk based on the sorption curve, considering the relative humidity of the air tending to $1$ and the operating conditions of the dryer. Nešić & Vodnik (1991) already warned about the reduction in the concentration of air saturation due to the presence of the solid phase in the drop. In practice, this would equate to an $X_{cr}$ greater than experimentally determined. Studies are underway for equating $X_{cr}$ as a function of the equilibrium properties of the solid-liquid-gas phases.

The droplet spatial temperature profile predicted by the model over time as a function of the normalized radius ($R_d\left(t\right)/R_{d,0}$) is shown in Figure 4. Initially ($t=0s)$, the particle is entirely at the initial temperature and, then, the heat rate received on its surface is transferred to its center gradually. This temperature tends to stabilize rapidly near the wet bulb temperature. This rapid heating can be seen in both Figures 3 and 4. The results obtained allow us to conclude that the temperature distribution inside the drop must be considered during drying, as suggested by Mezhericher et al. (2007). The particle shrinkage profile predicted by the model is illustrated in Figure 5, but indirectly it can also be seen in Figure 4. This shrinkage occurs predominantly at $23.0$°C for milk and $32.6$C for silica and during the rapid initial heating of the droplet, the rate of evaporated water is still reduced corresponding to the value of approximately $0.97$ of the initial droplet radius for both skim milk and colloidal silica.

Figure 4
Spatial temperature profile during drying in the radial direction: skim milk droplet (a); colloidal silica droplet (b).

Figure 5
Shrinkage profile of a single droplet of skim milk with drying air at 50C (a), of colloidal silica with drying air at 101C (b).

The perfect shrinkage result was already expected, as it was a consideration of the developed model (Equation 17), which is a coherent consideration for the first drying stage. However, the dimensionless Peclet number can be used to indicate the tendency of crust formation in a droplet and therefore can express the degree of deviation from the perfect shrinkage (Chew et al. 2015). Thus, the compaction of solids needs to be analyzed to obtain a more accurate transition criterion from the free water evaporation stage to the crust formation stage. For the first drying stage, shrinkage is proportional to water evaporation. That is, a good prediction of the droplet mass directly influences the droplet shrinkage mapping. As well as the temperature, RMSE and $R^2$ values are presented in Table IV for the mass variable prediction.

Table IV
Performance evaluation (RMSE [×10⁷kg] and R$R^2$) for droplet mass simulations.

Product	Correlation	Evaluated at T∞		Evaluated at $$\overline{T}$$
Skim milk	Ranz-Marshall	$RMSE=0.6216$	$R^2=0.9835$	$RMSE=0.7201$	$R^2=0.9779$
	Whitaker	$RMSE=1.1801$	$R^2=0.9407$	$RMSE=1.3388$	$R^2=0.9237$
Colloidal silica	Ranz-Marshall	$RMSE=2.4055$	$R^2=0.7762$	$RMSE=2.1906$	$R^2=0.8144$
	Whitaker	$RMSE=2.0124$	$R^2=0.8434$	$RMSE=1.8506$	$R^2=0.8676$

The simulation of skim milk drying obtained good agreement between simulated and experimental values for both correlations and reference temperatures with better results using Ranz-Marshall correlation. Especially for silica, using $T_{\infty }$ as a reference temperature, despite relatively acceptable $R^2$ values, the simulated models were greatly overestimating the water evaporation, which caused a large drop in the quasi-equilibrium temperature of water evaporation, resulting in low $R^2$ values for temperature (see Table III). The dynamic profiles of the droplet mass during the skim milk droplet drying (Figure 6a) and colloidal silica droplet drying (Figure 6b) presented were for the best conditions presented in Table IV, that is, the ones with the smallest errors, namely, Ranz-Marshall-$T_{\infty }$ for skim milk and Whitaker-$\bar{T}$ for colloidal silica.

Figure 6
Droplet mass time profile of single droplet drying of skim milk with drying air at 50C (a), of colloidal silica with drying air at 101C (b).

It can be seen in Figure 6 that there is agreement between the simulated droplet mass values and those obtained experimentally (deviation $<9\%$ for skim milk and $<7\%$ for colloidal silica). It is important to point out that in the experiments of Nešić & Vodnik (1991), the droplets were generated with a micro syringe with an initial mass variation of up to $3\%$. Droplets containing insoluble materials tend to have an evaporation behavior more similar to that of pure water, and after the first drying stage, the crust offers less resistance to water migration since there is no formation of agglomerates that significantly impact the mass transfer. For drops containing dissolved materials, with evaporation, there is an increase in the concentration of non-volatile components on the surface (Nešić & Vodnik 1991). However, for the first drying stage studied here, these differences were not observed in relation to the droplet components and the evaporation is mainly dependent on the air characteristics.

Evaluating the diffusivity of water in the particle during the spray drying process is another important factor that should be considered for better determination of $X_{cr}$. During drying, the solids diffuse to the center of the particle, as opposed to the diffusion of moisture. For milk, the lower diffusivity of fats and proteins leads to an increase in their concentrations on the sample surface, forming what is commonly referred as a film of milk cream (Kentish et al. 2005). This film formation adds a resistance to the diffusion of moisture through the surface due mainly to the fact that the fats have a hydrophobic characteristic. The effective diffusivity of water in skim milk was evaluated by Malafronte et al. (2015) and a mathematical correlation was suggested to describe the diffusivity as a function of water content and temperature. The effective diffusion coefficient is evaluated using a combination of nuclear magnetic resonance (NMR) and parameter estimation.

A strategy that can also be adopted to improve the performance of the proposed model is to seek to predict the variation of parameters through empirical correlations. Silveira et al. (2017) tested several parametric equations for the variable diffusion coefficient with the concentration, obtaining a significant improvement in the description of the process, when compared with works presented in the literature.

CONCLUSIONS

Temperature and mass history of the first stage of a single spherical droplet drying can be successfully predicted by the model presented. There were not significant differences observed for the first stage studied concerning to the droplet components (dissolved or insoluble materials) and the evaporation is mainly dependent on the air characteristics. The model was validated with experimental data for skim milk and colloidal silica drying, but it can be used for other materials.

The heating time of the droplet until reaching the quasi-equilibrium evaporation temperature is relatively short ($\Delta t\approx 7s$) for both materials and the evaporation during the first stage occurs mainly at the wet bulb temperature of the air ($T=23.0$°C for skim milk and $T=32.6$°C for colloidal silica). Discrepancy between simulated and experimental values did not exceed $9\%$ for skim milk and $7\%$ for colloidal silica during the first drying stage simulation. The best result was obtained for skim milk with the Ranz Marshall correlation evaluated at $T_{\infty }$ ($R^2=0.9835$). It is important to point out that the correlation is also a source of uncertainty for the simulations, as it deals with empirical parameter adjustments, which can reach $30\%$ of error in the correlations used, errors that are not so unusual for this type of correlation. The choice of the most suitable correlation and reference temperature can vary according to the characteristics of the air and the material to be dried. In this study, considering the applicability of the model in a more generic way, Whitaker correlation evaluated at the film temperature $\bar{T}$ showed better results, minimum of $R^2=0.7403$ for all simulations performed.

It is necessary to evaluate the value of $X_{cr}$, equating it according to the balance of the phases present, and the diffusivity of water in the particle during the entire drying process. A proposal to improve the performance of the proposed model is to predict the variation of parameters through empirical correlations, using, for example, a variable diffusion coefficient with concentration. Another alternative is to use the Peclet number as a transition criterion between the constant rate drying stage to the decreasing drying stage.

ACKNOWLEDGMENTS

The authors are grateful to CNPq – Conselho Nacional do Desenvolvimento Científico e Tecnológico – for the financial support.

• ABDULLAHI H, BURCHAM CL & VETTER T. 2020. A mechanistic model to predict droplet drying history and particle shell formation in multicomponent systems. Chem Eng Sci 224: 115713. doi: 10.1016/j.ces.2020.115713.
• BIRCHAL VS & PASSOS ML. 2005. Modeling and simulation of milk emulsion drying in spray dryers. Braz J Chem Eng 22: 293-302. doi: 10.1590/S0104-66322005000200018.
• BOTH E, SIEMONS I, BOOM R & SCHUTYSER M. 2019. The role of viscosity in morphology development during single droplet drying. Food Hydrocoll 94: 510-518. doi: 10.1016/j.foodhyd.2019.03.023.
• CAMPOS EC, CHAVES DHS, SOUZA PC, SOUSA MA, BIRCHAL MAS & BIRCHAL VS. 2021. Atomizers in spray dryers: a review. J Thin Films Coat Sci Technol Appl 8: 1-15. doi: 10.37591/jotcsta.v8i1.5517.
• CHAVES DHS, BIRCHAL VS, PASSOS ML & COSTA JR EF. 2019. Simulação da primeira etapa de secagem de uma partícula esférica considerando encolhimento e perfil de temperatura. In: Anais do XXXIX Congresso Brasileiro de Sistemas Particulados 2: 117064.
• CHEONG HW. 1983. The drying of small drops of particulate slurries. Doctoral’s thesis. Aston University, Birmingham.
• CHEW JH, WOO MW, CHEN XD & SELOMULYA C. 2015. Mapping the Shrinkage Behavior of Skim Milk Droplets During Convective Drying. Dry Technol 33: 1101-1113. doi: 10.1080/07373937.2014.985791.
• COSTA JR EF, VIEIRA RC, SECCHI AR & BISCAIA JR EC. 2003. Dynamic simulation of high-index models of batch distillation processes. Lat Am Appl Res 33: 155-160.
• DALMAZ N, OZBELGE HO, ERASLAN AN & ULUDAG Y. 2007. Heat and Mass Transfer Mechanisms in Drying of a Suspension Droplet: A New Computational Model. Dry Technol 25: 391-400. doi: 10.1080/07373930601184569.
• DONADEL K, FELISBERTO MD & LARANJEIRA MC. 2009. Preparation and characterization of hydroxyapatite-coated iron oxide particles by spray-drying technique. An Acad Bras Cienc 81: 179-186. doi: 10.1590/s0001-37652009000200004.
• FU N, WOO MW, SELOMULYA C & CHEN XD. 2013. Shrinkage behaviour of skim milk droplets during air drying. J Food Eng 116: 37-44. doi: 10.1016/j.jfoodeng.2012.11.005.
• GONZÁLEZ-TOXQUI C, GONZÁLEZ-ANGELES A, LÓPEZ-AVITIA R & MENDOZA-MUÑOZ I. 2020. Time and energy reduction on grape dehydration by applying dipping solution on freeze drying process. An Acad Bras Cienc 92: e20190072. doi: 10.1590/0001-3765202020190072.
• HANDSCOMB C, KRAFT M & BAYLY A. 2009. A new model for the drying of droplets containing suspended solids. Chem Eng Sci 64: 628-637. doi: 10.1016/j.ces.2008.04.051.
• KENTISH S, DAVIDSON M, HASSAN H & BLOORE C. 2005. Milk skin formation during drying. Chem Eng Sci 60: 635-646. doi: 10.1016/j.ces.2004.08.033.
• KESHANI S, DAUD WRW, NOUROUZI M, NAMVAR F & GHASEMI M. 2015. Spray drying: An overview on wall deposition, process and modeling. J Food Eng 146: 152-162. doi: 10.1016/j.jfoodeng.2014.09.004.
• LIN SXQ & CHEN XD. 2002. Improving the Glass-Filament Method for Accurate Measurement of Drying Kinetics of Liquid Droplets. Chem Eng Res Des 80: 401-410. doi: 10.1205/026387602317446443.
• MALAFRONTE L, AHRNÉ L, KAUNISTO E, INNINGS F & RASMUSON A. 2015. Estimation of the effective diffusion coefficient of water in skim milk during single-drop drying. J Food Eng 147: 111-119. doi: 10.1016/j.jfoodeng.2014.09.032.
• MALISKA CR. 2013. Transferência de calor e mecânica dos fluidos computacional. 2nd ed. Rio de Janeiro: LTC, 453 p.
• MEZHERICHER M, LEVY A & BORDE I. 2007. Theoretical Drying Model of Single Droplets Containing Insoluble or Dissolved Solids. Dry Technol 25: 1025-1032. doi: 10.1080/07373930701394902.
• NEŠIĆ S & VODNIK J. 1991. Kinetics of droplet evaporation. Chem Eng Sci 46: 527-537. doi: 10.1016/0009-2509(91)80013-o.
• RANZ WE & MARSHALL JR WR. 1952. Evaporation from drops Part I. Chem Eng Prog 48: 141-146.
• SANDLER SI. 2020. Chemical and Engineering Thermodynamics. 5th ed. Nashville: J Wiley & Sons, 1040 p.
• SCHUTYSER MAI, BOTH EM, SIEMONS I, VAESSEN EMJ & ZHANG L. 2018. Gaining insight on spray drying behavior of foods via single droplet drying analyses. Dry Technol 37: 525-534. doi: 10.1080/07373937.2018.1482908.
• SILVEIRA JC, DA COSTA AOS & DA COSTA JUNIOR EF. 2017. Modeling of basil essential oil extraction with variable diffusion coefficient. Eng Agrícola 37: 717-726. doi: 10.1590/1809-4430-eng.agric.v37n4p717-726/2017.
• VICENTE J, PINTO J, MENEZES J & GASPAR F. 2013. Fundamental analysis of particle formation in spray drying. Powder Technol 247: 1-7. doi: 10.1016/j.powtec.2013.06.038.
• WHITAKER S. 1972. Forced convection heat transfer correlations for flow in pipes, past flat plates, single cylinders, single spheres, and for flow in packed beds and tube bundles. AIChE J 18(2): 361-371. doi: 10.1002/aic.690180219.
• YILMAZ B, CAKMAK H & TAVMAN S. 2019. Ultrasonic pretreatment of carrot slices: Effects of sonication source on drying kinetics and product quality. An Acad Bras Cienc 91: e20180447. doi: 10.1590/0001-3765201920180447.

Publication Dates

• Publication in this collection
  14 Apr 2023
• Date of issue
  2023

History

• Received
  3 Mar 2022
• Accepted
  21 June 2022