Advanced Shrinking Particle Model for Fluid-Reactive Solid Systems

In the present work, the extended shrinking film model (ESFM) was applied to a reversible reaction in which a solid dissolves and reacts with a component present in the liquid phase. The model considers the reactive solid dissolution in the liquid phase and the diminishing of its radius with the reaction time. Furthermore, the liquid film surrounding the particle, through which the liquid component diffuses to react with the dissolved solid, is considered radius dependent; thus, the model is based on the mass balance equations derived for the solid surface, liquid bulk, and the liquid film. The model consists of two ODEs and a PDE, solved numerically with gPROMS ModelBuilder 4.0. It was demonstrated that the model can cover a wide range of operation conditions, and it shows a high degree of flexibility, allowing the application to several kinds of solid-fluid processes, such as esterification, gasification, and steam cleaning for the removal of dangerous and polluting gases (CO2, SO2) from the main process stream and NO capture.


INTRODUCTION
In the chemical industry, a lot of solid-fluid processes of high importance are applied, such as reduction of iron oxides, combustion of solid fuels, gasification of coal, and removal of pollutants. Moreover, more solid-fluid systems appear nowadays in chemical reaction engineering, due to the increasing use of biomass technology. The general approach for heterogeneous reactions, more specific for solid-fluid reactions, is based on the porosity of the solid. The reaction system, which considers a non-porous reactant, present fewer steps than a reaction system with a porous one. In the former case, the intraparticle diffusion is avoided, and the reaction system presents the diffusion of the component through the film because the external solid surface is present. When the chemical reaction rate is rapid, the diffusion is the rate-limiting step, so a one-step reaction system is present. Furthermore, the film thickness transfer between the phases is considered constant during the overall process as it is related to the mass and diffusion coefficients and calculated through the dimensionless Sherwood number (Sherwood et al., 1975). Customarily, the analysis of the reactions with negligible porosity, under isothermal conditions has been based on the shrinking core model (Levenspiel, 1972). Most kinetic models existing in literature to describe the solid-fluid reactions (Gibilaro et al., 1970;Williams et al., 1972;Park and Levenspiel, 1975;Jidaspow et al., 1976;Fan et al., 1977;Verbaan and Crundwell, 1986;Box and Prosser, 1989) are based on gas-solid heterogenous reactions. About the application of a reactive solid inside liquid media, few articles are present in literature with the application of the above-mentioned kinetic model (Lindman and Simonsson, 1979;Ranjan Jena et al., 2003). Moreover, some of them present simplifications on the kinetic model as negligible diffusion resistance toward the particle or the simple application of the film theory (Le Coenta et al., 2003). The present article is aimed to enhance the application of the film theory on solid-liquid reaction system with a non-porous solid reactant through the application of the extended film shrinking model (EFSM). The substantial difference compared to the simpler film theory is the film thickness dependence on the radial coordinate of the particle (Green and Perry, 2007), which shrinks during the reaction. The physical steps involved during the reaction of a solid immersed in a liquid are (1) partial dissolution of the solid component at the solid surface; (2) diffusion of the liquid component from the bulk toward the solid surface through the film surrounding the particle; and (3) chemical reaction in the liquid film and liquid bulk. Starting from the identification and consideration of the physical steps, the model was set up through mass balances for each phase involved in the reaction with the aim to applicate the present model to various reaction systems as demonstrated in our previous work Salmi et al., 2017).

MODEL DEVELOPMENT Extended Shrinking Film Model (ESFM)
The ESFM is based on the following hypotheses (Figure 1): • A non-porous solid reactant dissolves and reacts with the compounds presents in the liquid phase. • The reaction takes place simultaneously in the film surrounding the particle, and in the bulk phase, the faster the reaction is, the more dominant the film reaction becomes. • The film surrounding the particle is particle radius dependent (t 2 > t 1 > 0), indeed, the smaller the radius becomes as reaction proceeds, the thinner the film is (δ).
However, the film thickness depends not only on the particle radius but a contribution by the degree of mixing of the reaction system is considered, too. As it is well-known, the thickness of the stagnant liquid film around the solid particle diminishes as the stirring rate of the system is increased, leading to the decrease of the particle radius due to the consumption of the solid reactant, thus to a more dominant bulk-phase reaction. These conditions are achievable in laboratory-scale experiments, in which a high stirring efficiency is easily materialized. In the industrial scale, however, it is not always possible to reach diffusionfree conditions in the treatment of dissolving solid materials when liquids and film reactions become prevailing. During the reaction, the particle dissolves reducing its radius symmetrically, however, different solid particle shapes can be considered, such as spheres, cylinders, slabs, or realistic more irregular particle geometries. No reaction inside the particle occurs as either the solid is considered non-porous or the reactions sufficiently fast to avoid any diffusion resistance inside the particle. Furthermore, for the sake of simplicity, solid particles of equal size and isothermal conditions are considered. Evidently, a more rigorous approach could consider both particle size distribution and energy balance equations. The reaction system is simulated in a batch reactor with ideal mixing, thus allowing to apply uniform concentrations into the liquid bulk phase excluding local deviations in the system. To set up the model, the mass balance needs to be written into the different zones of the system in which the reaction takes place. As discussed above, after the solid dissolution the reaction occurs only into the liquid film surrounding the particle and into the bulk of the liquid phase, so three zones need to be accounted for the mass balance, these are the solid surface, the liquid film, and the bulk of the liquid phase. Dimensionless equations were derived to express the fundamental relations between the quantities appearing in the model and to perform the numerical computations in a rational way. In the model we have a couple of ordinary differential equations (ODEs), which are the dimensionless concentrations in the liquid bulk (17) and the dimensionless amount of the solid substance (23). Additionally, partial differential equations (PDEs) are included that represent the dimensionless profiles in the liquid film (5). The latter is solved by choosing a second-order centered finite difference formula for the spatial derivatives (the number of discretization points was 40) to convert the PDEs to ODEs. The model equations were solved numerically by using the advanced modeling tool gPROMS Model Builder v.4 software (Russo et al., 2015a,b). A little reminder about the derivatization of the equations used in the model are reported below; for more details, see our previous work .

Liquid Film: Mass Balances and Film Thickness
Considering an infinitesimal control volume inside the liquid film, between the solid surface and the liquid bulk, the mass balance for a dissolved component (i) can be written as: where N i denotes the molar flux, A is the surface through which the components diffuse, Γ i is the generation rate in the case of a single chemical reaction, and V is the volume of element. For more details, see the list of symbols.
The difference between the diffusion fluxes is (N i A), the amount of substances (n i ) is equal to c i V. The volume element V in Equation 1 is let to shrink ( V→0), and Equation 2 is obtained after straightforward mathematical substitutions: Equation 2 is presented in a very general form, where the shape factor (s) allows the application on different kinds of particle shapes. The shape factor s = 2 corresponds to a spherical particle, s = 1 a long cylinder, and s = 0 is a slab. For a real particle with a non-ideal surface, the shape factors can have non-integer values even exceeding 2 (Salmi et al., 2010). The non-ideality of the surface is typically assigned to surface aspects, such as cracks and craters as well as to limited porosity. Fick's law is applied to describe the molar flux (N i ), and after developing the derivative in Equation 2, the mass balance become: The radial coordinate changes between the value r = R on the solid surface to r = R + δ at the interphase between the liquid film and the liquid bulk. The generation rates (Γ i ) are expressed by a standard stoichiometric approach: where R denotes the reaction rate, ν i is the stoichiometric coefficient of component i in the reaction.
Dimensionless film coordinate [x = (r-R)/δ)], concentration (y i = c i /c 0 ) and time (θ = t/τ ) are introduced, where δ is the film thickness, c 0 the initial concentration, and τ an arbitrary reaction time. The mass balance equation for a component i in the film becomes: where D i τ /R 0 2 is the first dimensionless number defined (AD 1 ), representing the extent of the molecular diffusion compared to the initial particle radius.
The following boundary conditions are needed to solve the PDE: According the to the film theory, the film thickness is δ = D i /k Li . It is calculated by using the standard correlation relating the Sherwood number (Sh) to the Reynolds (Re) and Schmidt (Sc) numbers (Wakao, 1984). The application to a stirred tank allows an estimation of the Reynolds number with the Kolmogorov correlation (Temkin, 1977).
The following equation was derived for the film thickness (δ) : Some considerations based on the amount of substance present into the reaction system can adapt the equations into a form better suited to practical calculations. Considering a generic instant of time, the amount of substance is n j = n P ρ P V P /M P , where n P is the total number of particles in the system. Similarly, the initial amount of substance is n 0j = n P ρ P V 0P /M P . For a general geometry V P /V 0P = (R/R 0 ) s+1 . Finally, the relationship between the particle radius and the amount of substance is obtained, Inserting the relationship into the expression for film thickness in Equation 6, the following equation is obtained: Equation 8 reveals that the film thickness diminishes from the initial value to zero as the solid reactant (j) is consumed. The film thickness for a fixed chemical system can be diminished by stirring: with increasing the stirring speed, the energy dissipated per time (ε) increases and the film thickness decreases according to Equation 8. The ratio R/R 0 is obtained from Equation 7 A dimensionless film thickness is defined as follows: which can be compressed defining a dimensionless number (AD 2 = α) including both physical parameters and the dissipated energy of the stirring device:

Bulk-Phase Balances
The mass balance of the dissolved component (i) in the liquid bulk phase can be written as The liquid volume (V L ) is assumed to remain constant as declared above. After division by the liquid volume we obtain: where a P represents the total outer particle surface area-to-liquid volume ratio, and it is expressed as A P /V L . For a general geometry a P /a 0P = (R/R 0 ) s . The balance equation becomes now After inserting the relationship between the amount of substance and the radius in Equation 7, the balance equation becomes: The flux at the film-bulk interface is calculated from.
Finally, introducing dimensionless concentration and time, the liquid bulk balance equation is transformed to a dimensionless form: where D i a 0 τ /R 0 is the third dimensionless number (AD 3 ), meaning the ratio between the diffusivity of the component weighted by the specific surface area and the initial radius of the particle. As revealed, AD 1 and AD 3 are related to each other as function of similar quantities.
The flux at the film-bulk interface is calculated from

Solid-Phase Balances
The mass balance equation for a solid phase component (j) is given by where N j denotes the diffusion flux at the dimensionless coordinate x = 0, that is, at the outer surface of the particle. The total outer surface area is A = aV L = a 0 (R/R 0 ) s V L as discussed above. The Fick's law gives the flux, which in this case becomes The balance equation for a general particle geometry can be written as: After inserting the relationship between the amount of substance and the radius in Equation 7, the liquid bulk balance is obtained: The solid-phase balance equation is transformed to a dimensionless form by introducing the dimensionless time: where c 0 V L /n 0j is the fourth dimensionless number (AD 4 ).

Contribution Analysis and Modeling Approach
The contribution analysis of each phase allows to know in which extent the reaction is proceeds either in the film or in the bulk phase. The contribution of the film and bulk phases to the process is evaluated calculating the instantaneous contributions of each reactions in the film and in the bulk phase, respectively. The time to the shift between the film to the bulk-dominant reaction (t shift ) is measured as the inflection point of the relative contribution analysis of both film and bulk phases. The relative contributions of the film (F) and bulk (B) for the reaction can be calculated from the following equations: a more detailed derivation of Equations 24 and 25 is reported by Salmi et al. (2017). As a summary, the model to be implemented consists of Equations 5, 11, 17, and 23. The set of ordinary and partial differential equations constituting the model were solved numerically using gPROMS ModelBuilder v. 4.0 software, discretizing the film thickness with a second-order centered finite difference approach, adopting 50 points of discretization grid.

RESULTS AND DISCUSSION
In the present article the model was applied, first, to an elementary and reversible homogeneous reaction, in which a

B(s) → B Dissolution 2. A + B ↔ C Reaction
As the reaction occurs in both bulk and film phases, it is necessary to define the reaction rate expressions in each phase: For the reaction, a simple second-order reversible rate expression was adopted in which the term (c i ) and (c i ') are the concentration of the i-specie referred to the liquid film and the liquid bulk, respectively. The model was applied on several reaction conditions investigating the influence of either the operation conditions or the physical parameters on the behavior of the system. Table 1 shows the parameters used in each simulation; the shape factor (s) was set equal to 3 to simulate a realistic solid particle. In addition, to calculate a 0 , the molar mass of the solid was fixed to 0.1 kg/mol and the solid density to 2,000 kg/m 3 .
The characteristic time, needed to make our calculation time dimensionless, was chosen arbitrarily to be τ = 1s. Due to the high number of simulations (about 60) Table 2 shows the parameters ranges adopted in the numerical simulation tests. Figure 2 shows the trend of the dissolution time (t dis ) against the kinetics and equilibrium constant at different value of the first dimensionless number. The white zones indicate areas of the space of parameters not investigated. The absolute value of Frontiers in Chemical Engineering | www.frontiersin.org the dissolution time changes with the dimensionless number, passing from hundreds of seconds (AD 1 = 1 × 10 −3 ) to just 1 s (AD 1 = 1). In general, as AD 1 is a function of both D i and R 0 , an increase in AD 1 value can be justified by either a higher molecular diffusivity or a lower particle radius. In both cases, the diffusion resistances become less influential, leading to a decrease of the dissolution time, as revealed in Figure 2.
For the lower AD 1 (Figure 2A), the particle "lives" in the system for much time with respect to the other cases. The situation is more linear for Figure 2B in which the higher the reaction constant, the faster the dissolution time. In this case, the dissolution time is independent from the equilibrium constant; in fact, at each value of the latter, the dissolution time is the same.
About Figure 2C, the trend is in general analogous to the first two figures, the dissolution time decreases with the reaction constant; the higher the reaction constant the faster the dissolution time. The situation is different at low values of the equilibrium constant (<100 [m 3 /mol]); in fact, the dissolution time presents a maximum at reaction constant equal to 0.015 [m 3 /(mol s)], and it decreases both diminishing and increasing the reaction constant.
Finally, Figure 2D shows the situation with the higher dimensionless number (AD 1 = 1), in this situation the dissolution time is lower. The trend can be divided into two ranges. The first one is for an equilibrium constant higher than 10, and the dissolution time is faster, increasing either the reaction or equilibrium constant. Contrarily, at equilibrium constant values <10 [m 3 /mol], the dissolution time is faster either with increasing the reaction rate constant or decreasing the equilibrium constant. The slower dissolution time for this situation is at an equilibrium constant close to 10 [m 3 /mol] and a lower value of the reaction rate constant. Figure 3 shows the trend of the shift time (t shift ); this time suggests when the reaction pass from the liquid film to the liquid bulk or vice versa. As demonstrates in the previous paragraph, the contribution at the bulk or film is given by the instantaneous contributions of each reaction in each phase. The trend of the shift time is quite similar to the dissolution time displayed in Figure 2. The time requested for the shift from the film to the bulk, for the first dimensionless number equal to 1 × 10 −3 , is higher than in the other three situations in which the shift time is <1 s. In this situation (Figure 3A), the higher the reaction constant the higher the shift time; however, the fixed the reaction rate constant the shift time is quite similar increasing the equilibrium constant. Regarding Figure 3B, the situation is linear. The shift time increases with the reaction constant from 0.01 [s], at the lowest value of the former for every value of equilibrium constant, to 0.15 [s] at the highest value of the former for every value of the equilibrium constant. Figure 3C shows an appreciable change in the shift time only for the equilibrium constant value below 100 [m 3 /mol]. In fact, above this magnitude, independent of the reaction constant, the shift time is the lowest. On the contrarily, under 100 [m 3 /mol] the maximum shift time in localized at the reaction rate constant equal to 0.02 [m 3 /(mol s)] and moving either at lower or higher reaction rate constant the shift time diminishes. Finally, in Figure 3D, the situation is different. The shift time presents a maximum between two local minima present at both low equilibrium and reaction constant and high reaction constant and low equilibrium constant, respectively. the total contribution of the film to liquid bulk is due to the diminishing size of particle and consequently of the film thickness, as declared in the present model.
In Figure 4, the time to the shift was calculated parametric with the reaction rate constant. As revealed, the shift changes of two orders of magnitude passing from 0.6 [s] when k = 1 × 10 −3 [m 3 /(mol s)], to about 220 [s] when k = 1 × 10 −1 [m 3 /(mol s)]. Figure 5 shows the shift time with the diffusion of the icomponent through the film as a parameter; in this case, the time for passing from the film to the bulk is instantaneous in both situations. However, a higher order of magnitude is required for the slower diffusion rate. This result is expected; in fact, a higher diffusion allows a rapid contact between the FIGURE 7 | Dimensionless numbers (AD 1 blue symbols and line, AD 3 red symbols and lines) vs. the dissolution time (t dis ) at α = 0 (no stirring, empty symbol) and α = 10 (high stirring value, full symbol), for a non-reactive system. dissolved solid, and the compound in the liquid phase leads to a faster diminishing of the solid particle and film thickness. Consequently, the momentaneous contributions of the liquid bulk became higher than that of the film and the reaction pass from the latter to the former. Figure 6 represents the shift time with the equilibrium constant as parameter. In this case, the shift time has the same order of magnitude for all the simulations within the set. The difference between the equilibrium constant equal to 1 [m 3 /mol] and 10 [m 3 /mol] is very minor. Contrarily, the shift time for the equilibrium constant value equal to 0.1 [m 3 /mol] is higher. This means that at higher equilibrium constants, the system approaches an irreversible state; thus, the influence of the chemical equilibrium becomes negligible on the reaction rate, and the shift time moves to lower values as the transition between film to bulk phase becomes the fastest possible.
Finally, the model was tested on the system without reaction to evaluate the dissolution time (t dis ) of the solid particle with respect to the dimensionless numbers (AD 1 and AD 3 ). In Figure 7, AD 1 and AD 3 were plotted vs. t dis , and not vice versa, as the figure is more illustrative as such.
Recalling that AD 1 and AD 3 are both proportional to the molecular diffusivity and the reverse of the initial particle radius, the trends with the dissolution time are qualitatively similar. As revealed, by either increasing the particle radius or decreasing the molecular diffusivity, a corresponding increase of the dissolution times was calculated. This information is rather logical as bigger particles need more time to dissolve and react; the same occurs if the solute has a low molecular mobility within the stagnant liquid film. Moreover, the variations in AD 1 and AD 3 parameters look appreciable when increasing the dimensionless stirring parameter (α), even if the trends are logical, because by increasing the stirring rate, a lower dissolution time was estimated.

CONCLUSIONS
The ESFM was successfully applied on bimolecular solid-liquid reversible reaction kinetics. The simulations were based on a general liquid-solid reaction system in which a non-ideal shape factor of the solid was considered, revealing the flexibility and the power of the ESFM model. A detailed investigation was conducted studying the cases of reversible reaction and the analogous non-reactive system, demonstrating how the operation conditions and the physical parameters affect the performances of the system.
The reported investigation revealed that both the initial particle radius of the solid phase and the physical nature of the involved components, that is, the molecular diffusivity, are crucial factors to determine the time of shift between the regimes. Thus, grinding correctly, the solid phase could lead to faster kinetics. Moreover, the eventual presence of chemical equilibrium could reduce the overall rate of the process. The nonreactive system case confirmed the high influence of the initial particle radius and the molecular diffusivity on the dissolution time, a fact even more evident when the system is stirred at higher rates.
General hypotheses were applied to simplify the mathematical treatment; however, a more specific and appropriate model development could be considered to further improve the model, as introducing the particle size distribution, as well energy balances for highly exothermic or endothermic reactions.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/supplementary material.

AUTHOR CONTRIBUTIONS
VR: conceptualization, supervision, and formal analysis. HG: conceptualization and formal analysis. TC: software and writing-original draft. RT: software. TS: conceptualization and supervision. All authors contributed to the article and approved the submitted version.

ACKNOWLEDGMENTS
University of Naples Federico II is acknowledged for funding The international agreement between University of Naples Federico II (IT) and Åbo Akademi University (FI) that allowed the exchange between researchers of the two affiliations, making the realization of this work possible. The economic support from Akademi of Finland is gratefully acknowledged (Academy Professor's grant given to TS).