CPFD study of volatile liquid injection in a dense gas–solid fluidized bed

This study uses the Eulerian-Lagrangian multiphase particle-in-cell (MP-PIC) approach, adopted in the computational particle ﬂuid dynamics (CPFD) code, Barracuda Virtual Reactor (cid:1) , to simulate the injection of a gas–liquid spray into a gas–solid ﬂuidized bed. The model accounts for evaporation of volatile liquid in the form of individually moving droplets and thin ﬁlms distributed on the surface of the solid bed particles. Comparison with previously reported temperature measurements for varying spray nozzle gas-to-liquid ratios suggests that the model is capable of predicting the overall features of the gas temperature distributions within the injection region at injection velocities of less than 50 m/s. Locally observed discrepancies between experimental and simulated temperature proﬁles are likely related to the presence of wet particle agglomerates in the spray region, as suggested by a developed energy dissipation model that, combined with the MP-PIC model, predicts the location and extent of agglomeration in the ﬂuidized bed.


Introduction
Injection of a liquid into a fluidized bed and subsequent wetting and drying of bed particles is an essential process to a range of industries including oil and gas refining, food processing, pharmaceuticals manufacturing, and polymer production (McDougall et al., 2005;Boyce, 2018).Although the amount of liquid injected over time into the bed for these operations is small relative to the bed volume and mass, several concomitant effects of the wetted fluidization stage increase the complexity of an already complex fluidized bed behavior (Boyce, 2018).The complexity of the phenomenon of liquid injection has spurred a number of researchers to conduct experimental investigations to increase the comprehension of the liquid injection, dispersion, and evaporation phenomena within the bed.Ariyapadi et al. (2003), Bruhns (2002), and Cocco et al. (2013) all conducted experiments in dense bubbling fluidized beds to investigate the behavior of liquid jets.Although the experiments of Cocco et al. (2013) were conducted under non-evaporating conditions, Ariyapadi et al. (2003), Bruhns (2002), and Cocco et al. (2013) all found that atomization of the liquid feed to form a spray was greatly suppressed in the presence of the high solids loading of the bed.Hence, Bruhns (2002) concluded that for dense gas-solid fluidized beds, the liquid injection mecha-nism is based on a liquid jet cavity appearing downstream of the nozzle exit rather than a spray of single droplets interacting with the bed material.Bruhns (2002) further argued that the momentum of the liquid feed is the primary controlling factor in the liquid injection mechanism rather than the distribution of single droplets, as these are not present.
Additionally, they all concluded that agglomeration of wetted particles occurs within or near the boundary of the liquid jet, where the shear stress exerted on the bed material by the jet velocity is low compared to the core of the liquid jet.
The subject of agglomeration of wetted particles in dense beds has been experimentally studied by several authors (McDougall et al., 2005;Mohagheghi et al., 2014;Leclère et al., 2004).Mohaghegni et al. (2014) found that increasing the superficial gas velocity reduces the agglomerate formation and lifetime.McDougall et al. (2005) found that the agglomeration tendency is not only dependent on the superficial gas velocity of the bed, but also depends on particle sphericity and the properties of the liquid, specifically its surface tension and viscosity.In addition, Leclère et al. (2001) and Motlagh et al. (2017) both noticed that the particle porosity is a decisive factor for the agglomeration tendency, as bed particles with higher porosity are able to absorb more liquid, leaving less of it available at the surface for liquid bridging between particles.Leclère et al. (2001) conducted their experiments in a bed with a maximum temperature of 250 °C and Motlagh et al. (2017) carried out their experiments at ambient temperature.Hence, both sets of experiments were performed at low bed temperatures relative to the boiling point of the injected liquids such as water and ethanol.At elevated bed temperatures, it is likely that the evaporation rate will be higher relative to the transport rate of liquid into the pores of the bed particles, thereby reducing the significance of the porosity.
Although formation and breakup of agglomerates are important aspects of liquid injection into dense gas-solid fluidized beds, other aspects have also been investigated experimentally.Prociw et al. (2018) found that increasing the gas-to-liquid ratio of the spray nozzle, used to atomize the liquid feed, improved the liquid distribution within the bed.Apart from increased spray jet angle, increased penetration length, and increased solids entrainment into the liquid jet, Prociw et al. (2018) also argued that smaller droplets formed at higher gas-to-liquid ratios would enhance the liquid dispersion.Bruhns (2002) investigated the effect of the bed temperature on the evaporation of water and the temperature distribution in the vicinity of the nozzle.By increasing the bed temperature from 120 °C to 180 °C the drying rate of wetted solids was observed to increase.Nevertheless, the thermal jet penetration lengths, i.e. the distance along the injection line from the nozzle exit to the point where the temperature approximately equals that of the interior bed, for the different bed temperatures were found to be the same (Bruhns, 2002).Finally, the minimum fluidization velocity of the bed, and thus its potential to defluidize during operation, has been found to depend on volatility and physical properties (such as viscosity and surface tension) of the liquid as well as the Geldart group classification and porosity of the bed material (McLaughlin and Rhodes, 2001;Seville and Clift, 1984;Wright and Raper, 1998;Wright and Raper, 1999;Du et al., 2006).
Within the methodology of computational fluid dynamics (CFD), various modeling approaches have been proposed to handle the experimentally observed effects of injecting volatile liquids into dense gas-solid fluidized beds.Several authors (Li et al., 2012;Ahmadi Motlagh et al., 2019;Pougatch et al., 2012;Werther and Bruhns, 2004) have expanded the two-fluid Eulerian-Eulerian model formulation (Anderson and Jackson, 1967) to include an additional continuous liquid phase.In addition to modeling of the liquid evaporation, Li et al. (2012) and Motlagh et al. (2019) accounted for agglomeration, while Pougatch et al. (2012) Greek letters well as droplet breakup and coalescence.Werther and Bruhns (2004) argued that agglomeration effects can be neglected when porous bed materials are used, and rather accounted for the effect of reduced evaporation rate with particle porosity.These studies showed acceptable agreement with experimental data in terms of injection region temperature distributions (Pougatch et al., 2012;Werther and Bruhns, 2004), vapor concentrations (Werther and Bruhns, 2004), agglomerate sizes (Ahmadi Motlagh et al., 2019), and the dispersion of liquid throughout the bed (Pougatch et al., 2012).Compared to the Eulerian particle and liquid phase description, a Lagrangian approach has the distinct advantage of being able to account for size distributions of both droplets and particles.The latter is of particular interest to e.g.detailed agglomeration studies of granulators, which are often conducted using discrete element methods to include deterministic collision and agglomeration models (Boyce et al., 2017;Sutkar et al., 2016;Mikami et al., 1998;Goldschmidt et al., 2003;Heinrich et al., 2003).However, such models are highly time consuming and not feasible for CFD modeling of commercial-scale facilities.Less time consuming discrete parcel methods, either neglecting or stochastically modeling particle collisions, have been adopted to simulate injection, evaporation, and reaction of a liquid feed into dilute fast fluidized beds such as FCC risers (Behjat et al., 2011;Nayak et al., 2005).However, to the authors' knowledge, the work of Zhao et al. (2009) is the only publication where this modeling approach has been used to simulate volatile liquid spray injection into a dense gas-solid fluidized bed.Zhao et al. (2009) used the hybrid Eulerian-Lagrangian multiphase particle-in-cell (MP-PIC) model with submodels for collisional transfer of mass, momentum, and energy in wetted gassolid fluidized beds proposed by Snider (2001) andO'Rourke et al. (2009) to simulate the liquid injection experiments by Ariyapadi et al. (2003).Although the penetration length was used as the only means of validation against experimental data, the MP-PIC modeling scheme with a stochastic collision model used by Zhao et al. (2009) was found to yield reasonable agreements with experimental data.Thus, in contrast to the liquid injection mechanism for dense gas-solid fluidized beds proposed by Bruhns (2002), injection and subsequent evaporation of a spray of droplets appears to capture at least one relevant aspect of the physics in this case.Naturally, this raises the question of whether the MP-PIC model can be used to model volatile liquid injection into dense gas-solid fluidized beds over a wider range of operating conditions and capture trends of e.g.varying bed temperature and droplet size as reported in the literature.Hence, in the attempt to answer this question while adding to the scarce model validation available in literature, this study adopts the MP-PIC model to simulate the experiments conducted by Bruhns (2002).Simulated bed temperature profiles are compared with experimental ones under varying operating conditions, including liquid mass flow rate, gas mass flow rate, and test liquids to validate the model setup.Further investigations of the effect of varying operating conditions on the behavior of the liquid injection region and dispersion as well as the formation of agglomerates, evaluated using an energy balance model developed for the MP-PIC scheme, are carried out.

Governing equations
The gas-solid-liquid flow is solved using the computational particle fluid dynamics (CPFD) methodology, which includes an implementation of the multiphase particle-in-cell (MP-PIC) formulation of the particle momentum description (Snider et al., 2011).The approach uses a fixed Eulerian grid for solving the fluid mass, momentum, and energy equations and sub-grid Lagrangian pointmass particles for which exchange terms exist in the gas-phase governing equations to ensure mass, momentum, and energy balance.The liquid feed is treated as a spray of droplet parcels with assigned sizes and properties that are prone to change as they experience heat-up, evaporation, and contact with the bed particles.Consequently, the droplets are largely treated in a similar manner to the solid particles, thus leading to highly similar formulations of the mass, momentum, and energy conservation equations to those describing the solid particles, outlined by e.g.Snider (2001) or Pannala et al. (2011).The most significant difference in description of the droplets and the particles is rooted in the ability of the liquid to evaporate and distribute onto the surface of the solid particles.Assuming that the efflux of vapor from the liquid surface proceeds at quasi-steady-state conditions at sufficiently low velocities to be neglected in the momentum balance, focus can be dedicated to the formulation of the mass and energy continuum and particulate conservation equations, for which the liquid evaporation is most significant.
For constant solid particle mass and no production of gasses through chemical reactions, the gas-phase mass and energy equations are formulated according to Snider et al. (2011) as: and @ @t respectively.In Eqs. ( 1)-( 2) the subscripts g and l indicate the gas and liquid, respectively.hg ;q g ;h g ;u g ; and p are the gas volume fraction, density, enthalpy, velocity, and pressure, respectively.Viscous dissipation and energy sources from other than the particulate phase (droplet and particles) are neglected.qin Eq. ( 2) is the gas heat flux, defined as the product of the gas thermal conductivity and the gradient of the gas temperature, while _ q D is the enthalpy diffusion source term arising from molecules of different gases diffusing into each other and transporting both energy/heat and mass (Snider et al., 2011).
The two source terms d _ m l and S h in Eqs.(1) and ( 2) are of particular interest to the current work, as they account for the gas mass production rate during liquid evaporation and the conservative energy exchange from the dispersed phase to the continuous phase, respectively.The two terms will be dealt with separately in the following.

Mass balance
The vaporization of the liquid is treated similarly to chemical reactions, and any liquid that is evaporated is added to the gas phase (Snider et al., 2011).Remaining liquid on the particles may be transferred to nearby particles via collision terms, as discussed later.The gas mass source term arising from the change in mass of the liquid, for either a droplet or a wetted solid due to evaporation, is defined in terms of the particle distribution function f: The quantity f gives the probability per unit volume of having a certain number of particles with mass, velocity, and temperature within the multi-dimensional differential elementdm p du p dT p , at an instant of time.Thus, integrating the product of the probability distribution function f and evaporation rate over a region of mass, velocity, and temperature space gives the gas mass source contribution within this region in units of kg/(m 3 •s) due to evaporation from the liquid phase.
Within the framework of CPFD, liquid may exist as individual droplets or it may exist on the surface of a solid particle as a liquid film.Consequently, the subscript p can refer to either a pure droplet (with very small solid core), a wetted solid particle, or a dry particle with zero liquid content.The solid part of a numerical particle will be denoted by s while the liquid part will be denoted by l: The term f in Eq. ( 3) is inspired by the Boltzmann equation from the kinetic theory of gases, and describes the statistical behavior of a molecular system that is not in an equilibrium state (Bird et al., 2007).
The mass balance for the liquid in a droplet or a wetted particle is made up of two terms: one term accounting for loss of mass due to evaporation and one term accounting for loss or gain of mass due to collision with other particles: The reader is referred to the work of Zhao et al. (2009) for description of e.g.how the rate of change in liquid mass due to collisions ðdm l;coll =dt in Eq. ( 4)) is defined and numerically calculated in the CPFD framework.The evaporative mass transfer from the liquid to the gas phase due to convection is approximated by a gradient diffusion model: where dm l;evap =dt is rate of mass transfer due to evaporation,c m;l is the mass transfer coefficient of the liquid,M is the molecular weight of the vapor species, while c and c 1 is the molar concentration of the vapor at the liquid film surface and in the surrounding gas, respectively.For simplicity, only one vapor species is considered in the current formulations.The ideal gas law is used to relate the temperature and the vapor pressure, determined by the Antoine equation, to c; while c m;l is calculated from the Sherwood number using the Ranz-Marshall correlation (Snider et al., 2011;Ranz and Marshall, 1952).The surface area of the liquid A l available for heat transfer is expressed in terms of the total surface area of the spherical particle (6)-( 7): Hence, for a small solid particle relative to the total particle size, the liquid surface coverage partition coefficient b becomes small and most of the heat transfer to and from the gas phase occurs via the liquid film.Conversely, for a small amount of liquid the size of the solid particle is almost similar to that of the entire particle, yielding a value of b close to one and thus a small surface area of the liquid film available for convective heat transfer between the gas and the liquid.

Energy balance
As the liquid film is treated independently from the solid particle, individual energy balance equations are formulated for the liquid and the solid particle.The energy balance for the solid particle is written as (Snider et al., 2011): where c p;s is the specific heat of the solid particle,T s is the solid core temperature,T g is the gas temperature, k g is the gas thermal conductivity, and Nu gÀs is the Nusselt number for convective heat transfer between the gas and the solid particle modeled using the correlation reported by Fan and Zhu (1998).The exposed surface area of the solid particle A s is determined in accordance with Eq. (6) as A s ¼ bA p : The second term on the right hand side of Eq. ( 8) results from the conductive heat transfer between the solid particle and the liquid film.The direction of the heat transfer, i.e. the sign of the temporal rate of change of the solid core temperature, is dependent on the local (Eulerian) gas temperature, solid core temperature, and liquid film temperature.In the spray zone, heated solids will be the main source of heat to both the liquid and gas, while solids that have been cooled upon contact with the liquid will heat-up outside of the spray zone by means of the heated gas entering the system.The temperatures in the solid core and in the liquid film are assumed to be spatially invariant, i.e. the solid core and liquid film can at one point in time only possess one temperature each, which can be different from each other.Analogous to Eq. ( 8), an energy balance for the liquid, existing either as droplets or liquid film on wetted solids, is written as: The last term on the right hand side of Eq. ( 9) accounts for the heat used to evaporate the liquid film (Lefebvre and Mcdonell, 2017), where the heat of evaporation h vap is calculated using the Watson equation (Watson, 1943).Finally, the energy exchange between the particulate and the gaseous phase S h in Eq. ( 2) can be written as a sum of the source terms related to the liquid and solid: and The enthalpies h s and h l are calculated from the specific heat capacities and heat of formations (Snider et al., 2011).The interphase drag function D p is calculated using the Wen-Yu/Ergun equation proposed by Parker (2015).

Collisional transfer of mass, momentum, and energy
During liquid injection into the fluidized bed, the spray of droplets will penetrate the bed due to the high initial velocity and momentum.Wetting of solids occurs by means of collisions between the liquid droplets and the bed particles as well as between wetted bed particles themselves.Hence, transfer of mass, momentum, and energy arise from interparticle collisions, which are paramount in predicting liquid dispersion.In the original formulation of the MP-PIC model, particle collisions are not included in the description of the particle distribution function f [29].Particle contacts are only accounted for by means of a particle normal stress evaluated on the Eulerian grid and whose gradient results in particle acceleration that prevents particle volume fractions from exceeding the close-packed volume fraction.O'Rourke et al.
(2009) included a collision term on the right hand side of the transport equation for the particle distribution function: / j denotes independent variables of the distribution function exclusive of spatial location and time, i.e. masses, temperatures, velocities etc. of the liquid film and solid particle.Further description of the transport equation for f is outside the scope of this paper.A model similar to the BGK collision model of Bhatnagar, Gross, and Krook (BGK) (Bhatnagar et al., 1954) is adopted as a closure relation for the collision term of Eq. ( 13).The BGK model assumes that the effect of the molecular collisions force a nonequilibrium distribution function to return to a Boltzmann equilibrium distribution function and that the rate at which this occurs is proportional to the molecular collision frequency (Bird et al., 2007).In this case, the BGK model has been reformulated to fit within the MP-PIC framework, specifically the particle distribution function f, hence referred to as a BGK-like collision model.As the collisions in the liquid-solid fluidized bed ultimately infer an equilibrium distribution in particle velocities, temperatures, liquid distribution and composition as well as other particle quantities, the collision model of O'Rourke et al. ( 2009) specifies two BGK-like collision terms: where f eq;u is the distribution that results if particle velocities are equilibrated and f eq;t is the distribution function that results from full collisional equilibrium of all particle quantities; mass, temperature, velocity etc..The two collision time scales s u and s t are proportional to the time between particle collisions s coll (O'Rourke et al., 2009): which depends on local particle velocity fluctuations, the size of the particle, the close-packed volume fraction, and the local solids volume fraction (O'Rourke et al., 2009).The proportionality constants K u and K t in Eq. ( 15) are determined from two collision frequency constants K L andK S : K S is the collision frequency for pure solids w ¼ 0 ð Þand K L is the collision frequency for pure liquids w ¼ 1 ð Þ .It takes 1=K S and 1=K L collisions to achieve velocity equilibrium for solid particles and liquid droplets, respectively.Hence, the collision model can be adjusted by tweaking of these two constants.Zhao et al. (2009) used K S ¼ 0:05 À 0:10 and K L ¼ 0:10 À 0:20 for which the model predictions of the jet penetration length were found to be insensitive.In addition to Eq. ( 16), a decay ratio R is defined in order to calculate K u and K t (O' Rourke et al., 2009): By including the collision term in the transport equation for f, the probability distribution function is modified to account for particle-particle contacts, whether it is pure liquid droplets or wetted solids, thus allowing the model to statistically predict liquid film masses, temperatures, and compositions across the full range of particles in the system.The reader is referred to the work of O'Rourke et al. (2009) or Zhao et al. (2009) for more detailed descriptions of the collision term along with definitions of the equilibrium distribution functions.

Energy dissipation model for agglomeration prediction
As particle collisions are a required predecessor for wet agglomeration to occur, an energy dissipation model is proposed to determine the extent of wet agglomeration during the transient CPFD simulations of the liquid injection into a dense gas-solid fluidized bed.The solution methodology is based upon a balance between the kinetic energy of rebounding wet particles and the energy dissipated by the liquid bridge forces, as proposed by Ennis et al. (1991).A similar criterion for agglomeration has also been adopted by a number of researchers in recent publications on CFD simulations of wet agglomeration, see e.g.(Wu et al., 2018;Sun et al., 2020).
The modeling approach adopted in the present study assumes that local properties of the wetted particles can be averaged within a computational cell to reduce the complex problem to a binary collision between two wetted particles with a representative particle size, surface roughness, liquid film thickness, and velocity.The binary collision of two equally sized wetted particles is schematically outlined in Fig. 1.The uniformly wetted particles with equal liquid film thickness approach each other with constant speeds along the centerline.During the approach, the total separation distance 2h sp is reduced.The separation distance is measured from the smooth particle surface, excluding the average height of the roughness asperities h a : When the distance between the two particles is less or equal to twice the liquid film thickness, the films start to merge, and forces are exerted on the two particles.Upon contact between the solid cores (including h a ) of the wetted particles an inelastic collision occurs, where a fraction of the initial kinetic energy is lost during the head-on contact.The direction of the par- ticles is reverted on the centerline and a liquid bridge is formed between the particles.
In accordance with the model of Sun et al. (2020), the kinetic energy of the rebounding particle is calculated as: where e s is the dry coefficient of restitution accounting for the loss of energy during the elastic collision,m s is the solid particle mass and U p;rel is the average relative velocity magnitude of the particles.The liquid bridge formed during contact between the two films is responsible for dissipation of the kinetic energy within the system.The dissipative energy of the liquid bridge can be estimated by spatial integration over the liquid bridge force F lb from the point of collisional contact (average height of the surface asperities h a Þ to half the rupture distance h rup (Wu et al., 2018): where 2h sp is the varying separation distance between the two wetted particles.The critical distance between the particles during the rebounding beyond which the liquid bridge ruptures and no agglomeration of wetted particles takes places is calculated by (Wu et al., 2018;Sun et al., 2020): where h is the contact angle and V þ lb ¼ V lb =r 3 p is the dimensionless liquid bridge volume.Several approaches to estimating the liquid bridge volume have been used in literature (Wu et al., 2018;Sun et al., 2020;Darabi et al., 2010).For simplicity the expression of Wu et al. ( 2018) is adopted, as they argued that the separation distance has little effect on the liquid bridge filling during the collision, which can thus be neglected: In Eq. ( 21), d þ l ¼ d l =r s is the dimensionless liquid film thickness prior to collision.

Liquid bridge force
The liquid bridge force, required for solving Eq. ( 19), can be estimated as the sum of the dynamic capillary force and the static viscous force (Ennis et al., 1991;Pitois et al., 2000).Gravitational forces are neglected and only the normal force component of the liquid bridge force is included in the analysis, i.e. tangential and torsional forces are not accounted for (Boyce, 2018;Pitois et al., 2000).The model of Pitois et al. (2000) for the liquid bridge force is adopted: where X V is a correction factor introduced by Pitois et al. (2000) to reduce the viscous contribution to the liquid bridge force in the case of small liquid bridge volumes and larger separation distances: In order to determine the relative interparticle velocity U p;rel during the collision, the principles of kinetic theory of granular flow are adopted, for which the granular temperature, defined as ( (Zhao et al., 2009;Goldschmidt et al., 2002): relates to U p;rel in the following manner (Darabi et al., 2010): U p;rms in Eq. ( 24) is the fluctuating particle velocity, calculated as the mass average particle velocity deviation from the local (cell) mass averaged particle velocity ( (Zhao et al., 2009).For the sake of simplicity, the relative or collisional velocity is kept constant in Eq. ( 22).The contributions from the capillary and viscous force to the total force exerted by the liquid bridge on two approaching spherical particles at constant velocity and liquid volume is presented in Fig. 2a.Comparing the magnitude of E kin and W lb provides the basis for evaluating when wet agglomeration occurs, i.e. when W lb P E kin : By varying the particle wetness in terms of liquid-tosolid (L/S) mass ratio and the collisional particle velocity a critical agglomeration curve can be obtained, as plotted in Fig. 2b.For a given value of L/S andU p;rel , Fig. 2b enables evaluation of whether agglomeration occurs for specific liquid and particle properties.

Simulation setup and conditions
The bed temperature measurements reported by Bruhns (2002) were selected for validation of the liquid injection and evaporation model.Bruhns (2002) carried out experiments in a bubbling fluidized bed with a rectangular 1 m Â 0.5 m cross-section using both non-porous quartz sand and porous FCC catalyst particles as bed material, injecting either ethanol or water as test liquids.The experimental setup is depicted in Fig. 3. Electrical heating rods were placed above the gas distributor to heat the fluidized suspension to 153 °C.Experiments were carried out at ambient pressure.The superficial gas velocity was kept at 0.3 m/s.15 °C water or 20 °C ethanol was injected horizontally into the fluidized bed 1.1 m above the gas distributor (Werther and Bruhns, 2004).An externally mixed spray nozzle with a spray angle of 20°was used to atomize the liquid.Dimensions of the spray nozzle can be found in (Bruhns, 2002).
The CPFD implementation of the MP-PIC scheme in the commercial software Barracuda Virtual Reactor Ò 20.0.2 was used to simulate the process of liquid injection, dispersion, and evaporation in the three-dimensional bubbling bed of Bruhns (2002).Several geometric simplifications were made to the experimental setup of Bruhns (2002).The cyclone was neglected and elutriated solids were artificially recirculated back to the bottom of the bed to ensure constant solids inventory.The internal heat exchanger was disregarded and fluidization air of 153 °C was uniformly supplied to the bottom face of the computational domain, determined by the location of the gas distributor, confer Fig. 3. Non-porous Geldart B quartz sand with a Sauter mean diameter (SMD) of 122 lm was used as bed material.The bed height of 1.7 m reported by Werther and Bruhns (2004) for the FCC catalyst bed material was adopted for the quartz sand.
Liquid feed was injected as droplets with a predefined size range and velocity.The latter can be determined by means of momentum conservation when provided with the mass flow rate of the gas _ m g and liquid _ m l (Walzel, 1993): Eq. ( 26) assumes that the liquid is accelerated by the gas up to a weighted velocity magnitude of U inj in accordance with Walzel (1993).The expression is valid for subsonic conditions and neglects the momentum exchange with surrounding gas, and is thus expected to overpredict the actual injection velocity.Apart from correlations for the vapor pressure and the heat of vaporization, material properties of the injected liquids, specifically, viscosity, heat capacity, and thermal conductivity, were specified as polyno-mial functions with respect to temperature.The liquid density was kept constant and equal to the value at the injection temperature.Bruhns (2002) provided SMDs of the droplets generated in air at ambient conditions for the two-fluid nozzle at various gas-toliquid (GLR) flow ratios [m 3 /h/(L/min)].The droplet SMDs were reported to decrease with increasing GLR.The droplet size distributions, provided as boundary conditions for injection of liquid parcels, were estimated using a Rosin-Rammler distribution with a distribution parameter of q ¼ 2:5 applicable to most sprays (Lefebvre and Mcdonell, 2017).The particle size distribution of the bed material was defined in a similar manner with a minimum and a maximum particle size of 5 lm and 350 lm, respectively.Fig. 4 depicts the droplet size distributions for three GLRs.
Table 1 lists the main parameters related to the operating conditions used to validate the model.Simulations were carried out by initializing the bed material in the domain with a volume fraction of h p ¼ 0:5; slightly below the close-packed volume fraction h p;CP ¼ 0:58 to ensure a smooth start-up of the simulations.Transient simulations with a variable time-step were conducted by simulating the bubbling fluidized bed for 2 s without injection of liquid.The velocity of gas and liquid injected by the simulated nozzle was then linearly ramped from t ¼ 2s to t ¼ 3s; after which the velocity was kept constant.Air was used as atomizing gas for all cases at a flow rate specified in Table 1 determined under the conditions of injection.Simulations were conducted for 20 s with the last 10 s used for time averaging of quantities to ensure proper repeatability of the results.The time-step was automatically adjusted during the simulation to ensure a Courant number between 0.6 and 1.0.Using a Cartesian cut-cell, staggered grid arrangement totaling 488 000 cells with local refinements made to the injection zone, an average time-step ofDt ¼ 1:1Â 10 À4 -2:2 Â 10 À4 s, depending on the operating condition, was  obtained.The chosen mesh configuration yielded cell lengths of 8.5 mm, corresponding to 70 times the SMD of the bed material.4.2 million parcels were used to properly discretize the bed mass inventory.
K L and K S of Eq. ( 16) were set to 0.1 and 0.25, respectively, based on results of preliminary sensitivity studies.The ratio of the collision constants R provided in Eq. ( 17) was set to 1, per Barracuda Virtual Reactor Ò 20.0.2 default, which according to the sensitivity study also yielded the best fit with the experimental data of Bruhns (2002).Results of preliminary sensitivity studies, which also include effect of mesh size, parcel count, time-step, simulation time, droplet size, and injection velocity magnitude are provided in the supplementary material.

Model validation against temperature measurements
For validation purposes, experimental and simulated temperature profiles are compared for the operating conditions of Table 1.Fig. 5 and Fig. 6 depict the instantaneous gas temperature field for operating condition W-8.3-50 in the vicinity of the nozzle over the last second of simulation.Generally, the lowest temperatures are observed near the injection point on the wall, increasing from left to right towards the bed temperature of 153 °C.
Comparing the temperature fields on the horizontal (Fig. 5) and vertical (Fig. 6) plane through the injection point, a symmetric temperature distribution is observed on the horizontal plane, while the upward flowing fluidization gas bubbles tend to deflect the jet slightly upwards, as shown on the vertical plane.This behavior is also apparent from the distribution of droplets, confer Fig. 7 and Fig. 8.Most of the droplets lose their horizontal momentum relatively fast, after which they are carried upwards with the fluidization gas and dispersed into the bed by the mixing of the solid particles.Fig. 7 and Fig. 8 reveal that the droplet sizes are considerably reduced after a short period of time due to evaporation and distribution of the liquid onto the solid particles.Consequently, most of the droplets outside the main spray region have radii of less than 5 lm.
The spray jet undergoes an oscillatory behavior with time as passing bubbles allows the atomization gas and droplets to extend further into the bed.The expansion is followed by a subsequent contraction in the jet length as large gas bubbles are released from the spray cavity, simultaneously carrying away both droplets and solid particles from the near jet region, confer e.g.Fig. 8a-b.This results in a pulsating behavior of the jet, much in line with observations of e.g.Briens et al. (2010;2021).
Upon temporal averaging of the gas temperature field over the last 10 s of the simulation time, Fig. 9 is obtained.By time averaging the temperature data along the injection line (dotted line in Fig. 9), the measured temperature profiles of Bruhns (2002) can be directly compared to the simulated ones by the present CPFD model.
Fig. 10 compares experimental and simulated time averaged temperature profiles along the injection centerline.The simulated gas temperatures result from an energy balance between droplets, solid particles and liquid films, confer Eqs. ( 2), (8), and (9).Fig. 11 plots the cell averaged droplet, solid, liquid film, and gas temperature along the injection centerline.The cell averaged temperatures are evaluated every 0.2 s during the last 10 s of simulation (from t ¼ 10 s to t ¼ 20 sÞ and subsequently averaged to yield a single time averaged value.Consequently, less smooth profiles than the gas temperatures plotted in Fig. 10 are obtained in this manner.Droplets enter the bed at the specified injection temperature and Fig. 4. Cumulative size distribution function of liquid droplets as generated in air at ambient conditions with the externally mixed two-fluid spray nozzle at various gas-to-liquid (GLR) [m 3 /h/(L/min)] ratios.

Table 1
Simulated operating conditions of Bruhns (2002) along with key input parameters for the CPFD simulation setup.Operating case identifiers are provided according to the test liquid, followed by the gas and liquid flow rate, i.e. for operating condition W-4.2-50 water is injected with a gas flow rate of 4.2 m 3 /h and a liquid flow rate of 50 L/h, while for operating condition E-8.3-25 ethanol is injected with a gas flow rate of 8.3 m 3 /h and a liquid flow rate of 25 L/h.are heated up towards the boiling point of the water.In the process, liquid is transferred via collisions to the solids, which in turn lowers their temperature.The core temperature of the wetted solids and the liquid film is almost the same, suggesting fast heat transfer between the two.The liquid film temperature is observed to rise above the boiling point.Such unrealistically high liquid temperatures are numerical artefacts that arise when very thin liquid films evaporate completely within a time step.Hence, the amount of liquid film mass with temperatures above the boiling point is negligible.The low temperature droplets have a pronounced influence on the gas temperature in the vicinity of the injection point.
As the liquid evaporates, the less wetted and dry solids become dominant with the gas temperature increasing towards the bed temperature, as also observed in Fig. 10.Fig. 10a displays the effect of varying the gas mass flow rate of the nozzle, while keeping the liquid mass flow rate constant.The   SMD of the droplet size distribution is changed according to the sizes listed in Table 1.However, as observed from the results of the preliminary droplet size sensitivity study provided in Fig. S.10 in the supplementary material, little or no effect of changing the droplet size distribution was observed.Hence, changes in the temperature distributions under varying operating conditions are primarily attributed to the change in flow rates of the liquid and atomizing gas.From the experimental data of Bruhns (2002), increasing the gas mass flow rate by a factor of two when going from W-4.2-50 to W-8.3-50 slightly increases the gas-liquid   varying liquid flow rate of the nozzle at constant gas flow is compared for experimental and simulated data.Experimental data is adopted from Bruhns (2002).penetration length as observed by shifting the temperature profile to the right on the abscissa.In both cases, similar shapes of the temperature profiles are observed with an initial steep increase in temperature.In this region (x < 3 cm), the CPFD model predicts similar trend, as the temperature increases steeply from the injection temperature.However, slightly lower temperatures are predicted by the CPFD model for both operating condition W-4.2-50 and W-8.3-50.The measured temperature profiles reveal a plateau 3-6 cm downstream of the injection point, which according to Bruhns (2002) indicates the formation of wet agglomerates sticking to the thermocouple in addition to affecting the local liquid evaporation rate.Such effects are not considered in the CPFD model.Hence, the gradual increase in gas temperature with increasing downstream distance from the nozzle continues for the simulated results, gradually declining in rate of change.Downstream of the plateau (x = 6-8 cm), the experimental measurements reveal a steep increase in the rate of change of the temperature profile, followed by a gradual reduction in temperature rise toward the interior bed temperature 8-12 cm downstream of the injection nozzle.Such trend is not captured by the CPFD model, as evident from Fig. 10a.Additionally, the CPFD model predicts a very gradual and smooth transition from the jet temperatures to the bed temperature at x = 10-20 cm downstream of the injection point.Here, the experimental data suggests a much more clear-cut transition from the jet region toward the interior of the bed.
With respect to the initial rate of temperature increase upstream of the experimental temperature plateau and the thermal extent of the gas-liquid jet, acceptable agreement between the experimental and simulated data is found for operating condition W-4.2-50.However, the overprediction of the thermal extent of the gas-liquid jet by the CPFD model is pronounced for the higher gas flow rate at operating condition W-8.3-50.
Keeping the gas mass flow rate fixed at 8.3 m 3 /h and varying the liquid mass flow rate between 25 L/h (W-8.3-25) and 50 L/h (W-8.3-50)yields the temperature profiles of Fig. 10b.Qualitatively, the experimental temperature profiles of case W-8.3-25, W-8.3-50, and W-4.2-50 look similar with a steep increase in gas temperature immediately downstream of the nozzle exit, followed by a region of more or less constant temperature, before the temperature increases again toward the bed temperature.Bruhns (2002) noted that an increase in the liquid flow rate at a constant atomizing gas flow rate slightly decreases the penetration length of the temperature profile, which they attributed to the lower GLR resulting in a lower momentum of the gas-liquid jet.However, when comparing the temperature profiles of W-8.3-25 and W-8.3-50 in Fig. 10b, the difference in the extent of the near-constant temperature plateau is the most noticeable, being larger for W-8.3-25 with the higher injection velocity.
The CPFD model does not capture this behavior.The simulated temperature profiles reveal higher temperatures up to 11 cm downstream of the injection point for the W-8.3-25 case than the W-8.3-50 case, which suggest that the initial droplet momentum is of secondary importance compared to the lower liquid flow rate.Further downstream, the two simulated temperature profiles approach each other, which suggests that the liquid flow rate and thus the amount of liquid available for evaporation mostly affects the core of the spray region, while the nozzle gas flow rate determines the overall extent of the spray region.
In order to assess the influence of the properties of the injected liquid, the local gas temperature profiles are compared for injection of liquid water and liquid ethanol in Fig. 12a.The flow rate of the liquid was kept at 25 L/h with an atomizing gas flow rate of 8.3 m 3 /h, thus yielding GLR = 20 m 3 /h/(L/min).In the experimental data reported by Bruhns (2002), the penetration length of the temperature profile is reduced for injection of ethanol compared to water, arguably due to a lower heat of evaporation of ethanol causing it to evaporate faster.The agglomerates formed by liquid ethanol bridging tend to decay more easily and cannot penetrate the fluidized bed as well as the agglomerates formed from liquid water.This causes the difference in the extent of the plateau near the boiling point of the two liquids (up to 6 cm and 9 cm downstream of the injection point for ethanol and water, respectively).The lower temperature of the E-4.2-25 plateau is due to the lower boiling temperature of ethanol compared to water.
As the CPFD model does not account for effects of agglomeration, this trend is not captured numerically.However, higher temperatures are observed over the entire length of the spray region for operating condition E-8.3-25, which is in agreement with the behavior observed in the experimental temperature profiles downstream of the temperature plateaus.
Fig. 12b evaluates the effect of keeping the liquid flow rate of ethanol constant while varying the atomizing gas flow rate.In agreement with the measurements, the two simulated temperature profiles are similar and with a longer thermal penetration length observed for E-8.3-25, although significantly overpredicted.Higher gas temperatures are found for E-8.3-25 than for E-4.2-25 in the vicinity of the nozzle most likely due to higher atomizing gas and droplet momentum, distributing the evaporation over a longer distance while increasing the amount of hot solid particles being entrained into the jet, ultimately leading to higher temperatures in this region.Acceptable prediction of the overall thermal jet penetration length is obtained at low injection velocity for case E-4.2-25.
Based on the presented findings, the CPFD model is capable of capturing the thermal extent of the spray region for low GLR ratios around 5 m 3 /h/(L/min) and thus low jet velocities of less than approximately 50 m/s.For higher momentum jets, the low temperature region is convected too far downstream in the model, leading to an overprediction of the thermal jet length.The experimentally observed effects of increasing the nozzle gas-flow rate while keeping the liquid flow rate are qualitatively captured by the CPFD model.Effects of varying liquid flow rate as well as test liquid are not sufficiently well captured by the model to allow for acceptable validation.The reason for these discrepancies are most likely attributed to the fact that the CPFD model does not include effects of wet particle agglomeration, which appears to be decisive from the temperature distribution within the injection region, as depicted in e.g.Fig. 12a.The lack of a temperature plateau in the simulated temperature profiles might also be attributed to an underprediction in the evaporation rate of liquid from wetted solid particles, which would extend the region of evaporation, leading to a less pronounced temperature plateau.The validity of these hypotheses is investigated further in the following section.Finally, it should be stressed that since the thermocouple measurements are inherently affected by not only the local gas temperature, but also by collisions with solids and droplets, a direct comparison with simulated gas temperatures, as performed here, might be misleading.This underlines the complexity of the flow in question and the difficulty of carrying out suitable measurements for unambiguous validation of numerical model predictions.

Effect of operating conditions on liquid dispersion
Fig. 13a plots the time averaged mass of liquid contained in either pure droplets or on wetted solids for each of the simulated operating conditions.The liquid content in the bed is evaluated every 0.2 s during the last 10 s of simulation (from t ¼ 10 s to t ¼ 20 sÞ and subsequently averaged to yield a single time averaged value.Most of the liquid mass in the bed is present as liquid film on the solid bed particles in the order of 25-100 relative to the mass contained in pure droplets, depending on the operating condition.By increasing the GLR, the liquid mass in the bed is reduced with respect to both droplets and solids.This trend is observed using both water and ethanol as test liquid.Generally, less liquid is present in the bed when ethanol is used as test liquid, most likely due to the lower heat of vaporization and higher vapor pressure, causing the ethanol to evaporate faster. Little or no difference in the simulated temperature profiles of e.g.W-8.3-50 and W-8.3-25 was observed in Fig. 10b, despite clear differences in the amount of liquid present in the bed for the two operating conditions as seen from Fig. 13a.Hence, the importance of the liquid dispersion and subsequent evaporation for the temperature distribution appears to be limited.Fig. 13b evaluates the amount of liquid present in pure droplets and wetted solids along the injection line.Only wetted solids and droplets within a cross-section of 8:5mm Â 8:5mm(corresponding to the size of a computational cell) of the nozzle injection line are considered.Droplet and wetted solids are grouped in sections with lengths of approximately 2 cm along the injection line to smoothen the profiles.Peaks in the liquid mass contained in pure droplets are observed 2.5 cm downstream of the injection point, reaching zero at x ¼ 10 cm after a rapid decay.Closest to the nozzle exit, the mass of liquid contained in droplets exceeds that of the mass contained on solid particles.However, due to method of evaluation, this peak is not observed in Fig. 13b.The decay in droplet mass is due to Fig. 12.Time averaged gas temperatures along the horizontal distance from the nozzle.Effect of (a) varying test liquid and (b) varying atomizing gas flow rate during injection of liquid ethanol is compared for experimental and simulated data.Experimental data is adopted from Bruhns (2002).three mechanisms: (1) evaporation of liquid from the droplets, (2) transfer of liquid mass from the droplets to the solid particles via collisions, and (3) hydrodynamic mixing in the bed, which carries the droplets away from the injection line.The three mechanisms are not readily distinguishable from each other.However, Fig. 13b clearly shows that liquid is transferred to the solid particles in the vicinity of the nozzle where the droplet mass peaks or slightly downstream of this point.For operating conditions W-4.2-50 and E-4.2-25 the former is the case, while for the other operating conditions, peaks in the liquid mass contained on wetted solids is found at around x ¼ 5cm: For operating condition W-8.3-25 an almost constant level in the liquid mass on wetted solids is found for x ¼ 2:5 À 7cm: Hence, a high gas flow rate yields a more pronounced solids-free jet cavity where less liquid mass is transferred from pure droplets to solid particles close to the nozzle exit.Instead, the liquid transfer occurs further downstream and is likely to extend for a longer distance into the bed due to higher jet momentum dominating over the fluidization velocity.The higher injection velocities also promote faster evaporation of liquid from the droplets, reducing the magnitude of the peaks observed in Fig. 13b for higher GLR, in line with findings from Fig. 13a.The region of high mass contained on wetted solids corresponds well with the locations of the temperature plateaus in Fig. 10 and Fig. 12.This suggests that agglomerates sticking to the thermocouple in the experiments by Bruhns (2002) is likely to occur, which would reduce the measured temperature in this region to a value between the wet bulb temperature and the boiling temperature.The appearance of more durable wet particle agglomerates in this region would promote such behavior in the measured temperature profiles.As the simulated gas phase temperature is not affected by formed agglomerates, the model does not predict the temperature plateaus.
Having considered the overall distribution of liquid mass in droplets and wetted solids, attention is dedicated to the distribution of liquid on solid particles, as this process is of importance in e.g.catalytic cracking reactions, coking applications, and formation of wet particle agglomerates.Fig. 14 plots a time averaged cumulative liquid distribution among wetted bed particles for each operating condition.
The curves are obtained by calculating the ratio of liquid-tosolid (L/S) mass for each computational solid particle in the bed with a liquid fraction larger than zero.A cumulative sum of the L/S mass ratio is calculated and normalized by the maximum value.This is done for each time-step of evaluation (every 0.2 s) and time averaged over the last 10 s over simulation from t ¼ 10 s to t ¼ 20 s: A high L/S mass ratio indicates that the particle has a large film layer deposited on its surface, and vice-versa for a low L/S mass ratio.The narrower the distribution, the more uniform the distribution of volatile liquid among the bed particles.Again, the GLRparameter appears to be decisive, as Fig. 14 shows that the overall L/S mass ratio is reduced and its profile turned more uniform when GLR is increased, independent of whether the increase is achieved by varying the feed flow or atomization gas flow.This is the case for results obtained using either ethanol or water.This observation is in agreement with the findings of Prociw et al. (2018) that the liquid distribution within the bed is improved by increasing the GLR.Operating conditions W-8.3-50 and E-4.2-25 yield almost identical results, despite the obvious differences in flow rates, liquid properties, and amount of liquid mass present within the bed.The same applies to operating conditions W-8.3-25 and E-8.3-25, which are highly similar in terms of L/S mass ratio.
Naturally, the L/S mass ratio depends not only on the liquid mass contained in the bed but also on the number of wetted solid particles.Fig. 15 plots the cumulative wetted mass fraction of uniquely wetted bed particles, i.e. particles that have not been wetted previously, over the last 10 s of simulation.The wetted bed mass is normalized by the total bed mass inventory of 1057 kg to yield the mass fraction of the bed that has been wetted as function of time.For comparative purposes the wetted fraction of the bed is set to 0 wt% at t ¼ 10 s: At t ¼ 20 s; between 0.25 wt% and 3.5 wt% of the bed has at some point been in contact with the liquid, depending on the operating condition.At these low mass fractions, the amount of uniquely wetted particles increases at a nearconstant rate.Assuming that a perfect and continuous mixing of all particles in the fluidized bed occurs, the near-constant wetting rate is expected to continue until the entire bed inventory has at some point been wetted.However, in practice, the wetted rate of nonpreviously wetted particles is expected to decrease with longer simulation times as the wetted mass fraction of the bed approaches 100 wt%.
Fig. 15 reveals that the rate of bed wetting is increased by decreasing the GLR.This is the case when either water or ethanol is used as test liquid.However, when using ethanol as test liquid, significantly lower bed wetting rates are observed as the ethanol is prone to evaporate faster from the droplets before contacting with the bed particles.

Wet agglomeration prediction
Based on the previously obtained results it seems unlikely that the lack of plateaus in the simulated temperature profiles is due to an underestimation of the evaporation rate of liquid from wetted   particles.Rather, it seems more likely that the distinctive plateaus in the temperature profiles found in the experimental data by Bruhns (2002) are the result of wetted particles and wet particle agglomerates contacting the thermocouple probe as well as affecting the local liquid evaporation rate.In this respect, the performance of the agglomeration model presented in Section 2.4 can be evaluated.
Fig. 16 plots the distribution of solid particles on top of the relative particle velocity field at various instances in time during gasliquid injection for operating condition W-8.3-50.Blue-colored particles denote a high L/S mass ratio, while red particles are dry and non-wetted.Fig. 16 shows that solid particles are wetted downstream of the injection point and carried upwards with the fluidization gas and the produced water vapor.The injection of the liquid droplets and the gas creates a more or less particlefree cavity downstream of the injection point, which extends slightly beyond the region of high jet and particle relative velocities.The existence of a more or less solids-free jet cavity is in line with findings of Cocco et al. (2013).However, contrary to observations made by Cocco et al. (2013), solid particles are wetted both within the jet cavity due to entrainment and along its border.
From the instantaneous local L/S mass ratio and relative particle velocity, the agglomeration criterion is evaluated in accordance with Fig. 2b to yield Fig. 17, where cells that comply with the agglomeration criterion within five separate instances in time are marked in red.The value of the surface asperity h a ; required to calculate W lb in, is set to h þ a ¼ h a =r s ¼ 0:0077; based on an average roughness height of 0.85 lm for fine-grained quartz sand particles with a mean size of 0.22 mm, reported by Alshibli and Alsaleh (2004).The prediction of agglomeration does not affect the simulation itself, hence, two subsequent time-steps are independent of each other with respect to agglomeration.Fig. 17 compares the agglomeration tendency for operating condition W-8.3-50 and W-8.3-25.In both cases, the energy dissipation model predicts wet particle agglomeration.Hence, the hypothesis by Bruhns (2002) that the plateaus in the experimental temperature profiles plotted in Fig. 10 and Fig. 12 are the result of agglomerate formation is highly plausible.The agglomeration occurs horizontally downstream of the injection point, in the same locations as high L/S mass ratios are observed in Fig. 16.At the boundary of the jet cavity, high L/S mass ratios are obtained while the relative particle velocity, highly affected by the momentum of the gas-liquid jet, has reduced enough for sticking rather than rebounding of colliding wet particles to occur, thus forming an agglomerate.This result is supported by Bruhns (2002), Cocco et al. (2013), andAriyapadi et al. (2003), who all found agglomeration to occur downstream of the injection point and in the low-shear region of the jet, as Cocco et al. (2013) and Ariyapadi et al. (2003) formulated it.
This observation is true for both operating condition W-8.3-50 and W-8.3-25.However, the agglomeration tendency for operating condition W-8.3-25 is lower than that of W-8.3-50, as fewer cells are prone to agglomerate in the considered time frame.This is to be expected, as the liquid flow rate of W-8.3-25 is only 25 L/h compared to 50 L/h for operating condition W-8.3-50 at the same gas flow rate, which additionally results in a higher injection velocity for operating condition W-8.3-25, thus leading to less agglomeration, as observed.
As the plateaus observed in the experimental temperature profiles of Fig. 10 and Fig. 12, according to Bruhns (2002), are indicative of the formation of wet agglomerates, the location and extent of the plateaus can be compared to time averaged estimations of the locations of agglomerate formation in the simulated data.For each case, i.e.W-8.3-50 and W-8.3-25, the horizontal distance (x-direction) from the nozzle exit of all cells prone to agglomerate and lying on the injection line is evaluated every 0.2 s from t ¼ 10 s to t ¼ 20 s: Upon averaging of these values, an upper and a lower bound, defined in terms of the standard deviation related to the average x-location at which wet agglomeration is predicted to occur, is calculated.The calculated ranges are presented in Fig. 18 alongside the experimental temperature profiles for the two investigated operating conditions.
The predicted agglomeration regions are compared to the location and extent of the temperature plateaus, which is the only experimental basis for comparing the predicted locations of agglomeration to what occurs in the fluidized bed.It is expected that agglomerates are formed in the first part of the plateaus and convected downstream with the atomizing gas, but also that agglomerates will form across the entire extent of the plateau.For operating condition W-8.3-50, the agglomerates along the injection line form approximately 10 cm downstream of the nozzle, while the agglomerates form closer to the nozzle (approximately 7.5 cm from the nozzle) for operating condition W-8.3-25.As previously noted, the measured temperature profiles suggest that the plateau is extended and shifted further downstream due to the increased injection velocity for operating condition W-8.3-25 compared to operating condition W-8.3-50.This behavior is not observed in the model predictions, where the agglomeration regions for the two tested operating conditions are of similar extent.However, the predicted locations of agglomeration are indeed located close to the expected ones and for operating condition W-8.3-25 the region of the predicted agglomerate formation even overlaps with the expected region deduced from the mea- surements by Bruhns (2002).This then provides some confidence in the model predictions, keeping in mind the large uncertainty related to validating the agglomeration method in this manner.
To quantify the degree to which agglomeration occurs, the maximum stable agglomerate size d agg is adopted from Darabi et al. (2010): where U Ã p;rel is the critical relative particle velocity derived from the critical agglomeration curve (confer Fig. 2b) at a specific solid par-ticle diameter d s ; obtained from interpolation of particle sizes to the cells.Predicted agglomerates with radii larger than the length of the local cell are disregarded to filter out unrealistically large agglomerates.Calculating the maximum stable agglomerate diameter for each of the cells prone to experience agglomeration at each evaluation time-step and summing their mass contributions within six sieve sizes yields Fig. 19.Hence, Fig. 19 is obtained by assuming that formed agglomerates do not breakup and is thus expected to overpredict the mass of solid particles contained in agglomerates after 10 s of evaluation.As expected, more of the bed material is predicted to be captured in agglomerates for operating condition   W-8.3-50 compared to W-8.3-25, suggesting that the simple way of predicting the agglomerate size has merits.For both operating condition W-8.3-50 and W-8.3-25 most of the total mass of particles captured in agglomerates is located in large agglomerates.This trend is in opposition to findings of e.g.Motlagh et al. (2019), where most of the mass is contained in smaller but more frequent agglomerates.The difference in the two model formulations lies in the feedback of the model, which in the current formulation is nonexisting, thus lacking the ability to account for agglomerate breakups and effect on mass and heat transfer limitations that affect the agglomerate diameter distributions.

Conclusions
Simulations of horizontal gas-liquid injection into a dense gassolid fluidized bed were carried out using the Eulerian-Lagrangian MP-PIC formulation implemented in Barracuda Virtual Reactor Ò .The computational particle fluid dynamics setup allowed for injection of a spray of liquid droplets into the bed, subsequently undergoing evaporation and distribution onto the hot bed particles by means of collisions.Previously reported experimental temperature data obtained using both water and ethanol as test liquid in a bed of non-porous quartz sand (Bruhns, 2002) was used as basis for validating the model under varying spray nozzle flow rates of gas and liquid.Good agreement was found at low injection velocities, i.e. high nozzle liquid mass flow rates and low atomizing gas mass flow rates.At high gas mass flow rates the extent of the thermal jet penetration length was overpredicted, independent of the test liquid.Experimentally observed variations in temperature profiles with test liquids and liquid flow rates were only observed to some extent in the simulated data, most likely due to the Lagrangian modeling of individually moving droplets and the lack of wet particle agglomeration effects.The model performed well in predicting the transfer of liquid from pure droplets to solid bed particles, as a larger nozzle gas-to-liquid mass ratio was found to provide a more uniform distribution of the injected liquid on the bed material, which is of importance to catalytic cracking reactions and wet particle agglomeration.
An agglomeration prediction model was presented, which is unprecedented within the framework of MP-PIC modeling.In agreement with expectations, agglomeration was predicted to occur under the simulated conditions with most of the agglomerate formation occurring in the low-shear regions of the gas-liquid jet where the spray cavity contacts with the bed material and where most of the particles are wetted.Less degree of agglomeration was satisfactorily predicted for a lower nozzle liquid flow rate and higher injection velocity using both qualitative and quantitative measures, the latter in terms of a maximum stable agglomerate diameter.
As agglomeration appears to be of importance in predicting not only the distribution of liquid throughout the bed but also the extent and behavior of the gas-liquid injection zone, further research is required to expand the current agglomeration model through integration with Barracuda Virtual Reactor.The model should entail a proper feedback loop to enable predicted agglomerates to undergo both growth and breakup, while affecting bed hydrodynamics as well as mass and heat transfer, which is not possible in the current post-processing format of the model.

Fig. 1 .
Fig. 1.Schematic representation of two colliding wetted particles of equal size and liquid content.Inspired by Ennis et al. (1991).

Fig. 2 .
Fig. 2. (a) Capillary and viscous forces of the liquid water bridge forces acting on two approaching 122 lm spherical and wetted particles obtained by solution of equation (22) for a constant relative velocity of 0.1 m/s.(b) Critical agglomeration curve derived from same conditions across a range of liquid-to-solid (L/S) mass ratios and relative particle velocities by setting equation (18) equal to equation (19).

Fig. 3 .
Fig. 3. Schematic of experimental fluidized bed setup used by Bruhns (2002) and the implementation of the setup in Barracuda Virtual Reactor.

Fig. 5 .
Fig. 5. Snapshots of instantaneous gas temperature distribution during gas-liquid injection for operating condition W-8.3-50.Sequence of contour plots on a horizontal x-y plane through the injection point is shown in time steps of 0.2 s from (a) to (e).Color bar is truncated at 85 °C and 150 °C.

Fig. 6 .
Fig. 6.Snapshots of instantaneous gas temperature distribution during gas-liquid injection for operating condition W-8.3-50.Sequence of contour plots on a vertical x-z plane through the injection point is shown in time steps of 0.2 s from (a) to (e).Color bar is truncated at 85 °C and 150 °C.

Fig. 7 .
Fig. 7. Snapshots of instantaneous location of droplets and solids (marked in yellow) during gas-liquid injection for operating condition W-8.3-50.Droplets are colored by their radius.Droplet and solids particles within a thin horizontal x-y slit centered around the injection point is shown in time steps of 0.2 s from (a) to (e).Color bar is truncated at 5 °C and 100 °C.

Fig. 8 .
Fig. 8. Snapshots of instantaneous location of droplets (marked in blue) and solids (marked in yellow) during gas-liquid injection for operating condition W-8.3-50.Droplets are colored by their radius.Droplet and solids particles within a thin vertical x-z slit centered around the injection point is shown in time steps of 0.2 s from (a) to (e).Color bar is truncated at 5 °C and 100 °C.

Fig. 9 .
Fig. 9. W-8.3-50 time averaged temperature distributions on (a) horizontal x-y plane and (b) vertical x-z plane through the injection point.Based on 10 s of averaging.Dotted line indicates the injection line on which time averaged data is extracted to yield temperature profiles.Color bar is truncated at 85 °C and 150 °C.

Fig. 10 .
Fig.10.Time averaged gas temperatures along the horizontal distance from the nozzle.Effect of (a) varying atomizing gas flow rate at constant liquid flow rate and (b) varying liquid flow rate of the nozzle at constant gas flow is compared for experimental and simulated data.Experimental data is adopted fromBruhns (2002).

Fig. 11 .
Fig. 11.Simulated cell and time averaged temperature distributions along the injection centerline.Obtained for case W-4.2-50.

Fig. 13 .
Fig. 13.Time averaged distributions of liquid mass in droplets and on wetted solids for different operating conditions.(a) depicts the total liquid mass appearing in the bed on either wetted solid particles as film (solid colored bars) and in pure liquid droplets (hatched bars) at an instant in time.(b) plots the sum of liquid present in droplets or wetted particles within the computational cells located on the injection line as function of the distance from the nozzle exit.

Fig. 14 .
Fig.14.Time averaged cumulative liquid distribution plotted against the liquid-tosolid (L/S) mass ratio for all wetted solid bed particles in the fluidized bed.

Fig. 15 .
Fig. 15.Cumulative mass fraction of bed particles that at one point have been wetted by liquid plotted against simulation time.The plotted mass fraction is normalized by the bed mass inventory (1057 kg) and set to 0 wt% at t ¼ 10 s for comparative purposes.

Fig. 16 .
Fig. 16.Snapshots of instantaneous location of all wetted solid particles during gas-liquid injection for operating condition W-8.3-50.Solids are colored by their L/S mass ratio while snapshots of instantaneous relative particle velocity is plotted on a vertical center plane colored by black in regions of high velocity (truncated at 10 m/s).Particles within a thin vertical x-z slit centered around the injection point is shown in time steps of 0.2 s from (a) to (e) during the last second of simulation.

Fig. 17 .
Fig. 17.Prediction of computational cells prone to agglomerate (red) during gas liquid injection for (a)-(e) operating conditions W-8.3-50 and (f)-(j) W-8.3-25.Sequence plots on a vertical x-z plane through the injection point is shown in time steps of 0.2 s from (a/f) to (e/j).Point of injection is marked by the black arrow.

Fig. 18 .
Fig. 18.CPFD estimated one-dimensional regions (along nozzle injection line) of wet agglomeration formation compared to mean temperature profiles along the distance x from the nozzle injection, measured by Bruhns (2002).

Fig. 19 .
Fig.19.Comparison of distribution of mass captured in agglomerates with respect to the predicted maximum stable agglomerate diameter.
included submodels to account for variation of droplet diameter within the continuous description of the liquid phase as