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.


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 (D/Dx) involves the specific mass  and specific enthalpy h with velocity u r (considering only radial velocity), which would encompass mass transfer within the particle (Equation 2).
There is no internal energy generation ( 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 ( is also zero).To express heat transfer, div ⃗ 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 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 ).

𝜕 (𝜌
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.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  and , 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.
The Equation 8describes the energy balance in each internal volume of the droplet.H represents the enthalpy of the control volume and ṁ the mass transfer rate.The sub-index w is for water and d is for droplet, 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).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,  -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.
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 16describes that heat enters the droplet by convection and mass evaporates from the outermost spherical layer of the droplet to the gaseous medium being  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 ṁ w between each of the various internal control volumes is proportional to the mass rate of water evaporated on the surface ṁv and to each volume in relation to the droplet total volume.
It should be noted that the last two terms of the Equation 16correspond 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(Δr 2 ), where Δ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 Δ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.
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 ∞ noting that the correlation is valid for Reynolds numbers between 0 and 200.
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 × 10 4 .All properties are evaluated at the air temperature T ∞ except  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.
The specific heat of the droplet considering it as an ideal mixture of solids and water can be evaluated by Equation 21where 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.
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  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  d is calculated as a function of moisture content X (Equation 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. (28) The values of water vapor densities on the surface of the droplet,  v,surf , and in the surrounding fluid,  v,∞ , 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.
v,∞ = H air P atm M air R (T ∞ + 273.15) (30) 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.
The correlations used for air properties are presented in Equations 34 to 37. The correlations for  air ,  air and k air were taken from Cheong (1983) 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 (SiO 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 SiO 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  d,w , k w , c p,w and  remain the same as those used for the skim milk drying simulation (see Table I).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.

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 ∞ ) 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 T (also called film temperature) between the drying air temperature T ∞ 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.
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.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 ∞ for skim milk and Whitaker-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 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 T indirectly included the influence of the boundary layer.
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 ≈ 150s for skim milk and in t ≈ 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 10g/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 (t)/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.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.×10 7 ×10 7 kg] and R 2 R 2 R 2 ) for droplet mass simulations.

Product
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 ∞ 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 ∞ for skim milk and Whitaker-T for colloidal silica.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.

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

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

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

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

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

Figure 6 .
Figure 6.Droplet mass time profile of single droplet drying of skim milk with drying air at 50°C (a), of colloidal silica with drying air at 101°C (b).
and c p,air from Sandler (2020).