Continuum Modeling of Slightly Wet Fluidization with Electrical Capacitance Tomograph Validation

: Gas–solid fluidized bed reactors are widely used in the power generation industry. The critical effect of the presence of liquid phase, either as a result of heat, chemical reaction or physical interaction, on the hydrodynamics of the reactor is well recognized by academic researchers and industrial operators. However, theory and simulation frameworks to predict such a condition using the continuum modeling approach are not yet available. This study first shows the significant changes in the flow pattern and distinguishable flow regimes in a slightly wet fluidized bed recorded by an advanced imaging technique. The study then describes the development and implementation of new mathematical formulations for wet particle-particle interactions in a continuum model based on the classic kinetic theory of granular flow (KTGF). Quantitative validation, carried out by comparing the predicted and measured fluidization index (FI) expressed in terms of pressure drop, has shown a good match. The prediction also demonstrates increased bubble splitting, gas channeling, slug-ging fluidization, and energy dissipation induced by liquid bridges developing from wet particle interactions. These characteristics are similar to those commonly observed in the fluidization of cohesive powders. This model constitutes an important step in extending the continuum theories of dry flow to wet particle-particle interactions. This will allow accurate description and simulation of the fluidized bed in its widest application including power generation systems that involve wet particle fluidization.


Introduction
Processing and handling of wet solids arise in a wide range of industrial processes, including chemicals, energy, pharmaceuticals, food, and minerals.For example, in power generation, oil/tar is formed in low-temperature regions of biomass thermal conversion reactors.Similarly, liquid may be present in fluid coking of tar sand (also called bituminous sand) or as a result of melting during particulate processing.In physical processing, liquid is commonly sprayed into mixers or fluidized beds for particle coating and granulation.The liquid presence in solid processing systems results in changing the inter-particle forces and energy dissipation at collisional contacts.The classical continuum model of granular flow, with constitutive relations adopted from the kinetic theory of granular flow (KTGF), has been widely used for many years to predict gas-solid flow behavior in dry conditions [1].This is commonly referred to in the literature as Eulerian-Eulerian or two-fluid modeling approach.The basic idea in the classic KTGF is that the particles are homogeneously distributed within an ensemble and can be defined as a continuous phase with stresses arising from random fluctuations, collisions and frictions, with each assumed to be dominant at a given flow regime.In analogy with thermal temperature in dense gases, a "granular temperature" Energies 2024, 17, 2656 2 of 20 is assumed to arise from the particle velocity fluctuations.The energy dissipation during collision is approximated by a restitution coefficient, taking into consideration the particle inelasticity.The random fluctuations and collisional stresses are approximated using equations derived from the Boltzmann equation while the fractional contacts are obtained from the Mohr-Coulomb law.The latter came as part of the efforts to extend the classic KTGF to slow frictional flow and, hence, constituted an important step towards taking this modeling approach to wider applications.In the case of a slightly wet condition, the particle-particle interactions are much more complex due to the formation of a liquid layer at the particle surface, which subsequently results in changing the overall hydrodynamics.
Geldart [2,3] classified the particles into four main Groups, A, B, C and D, according to their flow hydrodynamics when fluidized in air.The distinct characteristic behavior of each particle group was mainly attributed to the differences in their inter-particle forces.Following this, several studies have shown that the addition of a liquid to a dry bubbling bed may shift the hydrodynamics towards the distinctive behavior of cohesive particles (e.g., from Geldart Group B to A). Seville and Clift [4] and McLaughlin and Rhodes [5] experimentally studied the effect of the presence of small quantities of liquid (non-volatile silicone oil) on the fluidization behavior of glass beads belonging to Group B. Both studies have confirmed a shift in the hydrodynamic behavior towards another particle group; this shift was shown to coincide with the change in inter-particle forces.However, the change in bubble characteristics, the main factor that determines the bed hydrodynamics, was never discussed.The experimental work conducted by Makkawi [6] and Makkawi et al. [7] on a slightly wet bubbling fluidized bed of glass beads indicated critical changes in the bubble characteristics and clear channeling and voidage-wave instability as a result of liquid addition.Tardos and Pfeffer [8] reported an experimental study in which a fluidizing gas was found to chemically react with the solids producing new softer components and melting.The particle cohesion and stickiness were observed and bubble characteristics (velocity, size, and frequency) and minimum fluidization velocity were all found to change significantly as the reaction progressed.
One of the earliest attempts at the simulation of wet particles in a bubbling fluidized bed was reported by Mikami et al. [9] using discrete element modeling (DEM) coupled with computational fluid dynamics (CFD-DEM).The liquid bridge forces at the single particle level were analyzed based on the bridge shape and volume at the interface.The same approach was then followed by several researchers [10][11][12][13][14] and these are reviewed as follows.Xu et al. [10] used the CFD-DEM coupling model to study the effect of liquidinduced forces on gas-solid flow.It was found that particles tend to form agglomerates, resulting in an enlarged stagnant zone and increased fluctuation of bed pressure drop.Jian et al. [11] introduced a discrete model for gas-solid flow with wet (cohesive) granules.The transition from free-flowing to cohesive behavior was found to occur at a distinct value of the granular capillary number.Lim et al. [12] used the CFD-DEM coupling to study the mixing behaviors of wet granular materials in gas-fluidized bed systems.It was noticed that wet particles have a high tendency to form large aggregates with limited motion of individual particles.He et al. [13] used DEM to study wet cohesive particles in the presence of liquid bridges in a gas-fluidized bed.Symmetrical profiles of the velocity and granular temperature were noticed at low liquid contents, which then shifted to asymmetrical at a high liquid content.Most recently, Wang et al. [14] reported a CFD-DEM model to study the flow characteristics of dry and wet particles in a circulating fluidized bed riser.It was reported that the wet particles move in a single particle as well as an agglomerate.Most recently, Motlagh et al. [15] presented a comprehensive Eulerian-Eulerian CFD model to simulate evaporative spray injection into a hot gas-solid fluidized bed.The mode treated the three phases (gas, particle, and liquid) as Eulerian phases with a population balance model added to account for the particle size change after granulation.To account for the possibility of agglomeration, an energy balance model was used to assess the outcome of collisions between wet particles.This model is probably the first to treat wet fluidization in such a continuum modeling approach.In a more fundamental approach, Berzi1 and Fraccarollo [16] defined new closure equations to describe the fluid and particle stresses in a collisional suspension of particles and water flow.The model was strictly assumed to be valid for a situation where the particle inertia is dominant.Another major assumption was that the particle pressure scales well with the granular temperature, hence, the use of kinetic theory also permits to modeling of the granular viscosity.Accordingly, the effective fluid viscosity was assumed to arise from the combined effect of turbulence and a granular-like viscosity.The latter was modeled by simply replacing the particle density in the KTGF expression with a modified density taking into account the added fluid mass.The study has shown the feasibility of extending the two-phase continuous mathematical models beyond the boundaries imposed by the classic KTGF, which is in agreement with the outcome of the approach adopted in this study.

Motivation and Objectives of the Study
Although many of the industrial gas-solid flow systems are operated at slightly wet conditions, as described in the introduction, reliable theoretical models to predict the effect of the liquid presence in practical industrial-scale fluidization are currently not available.There are many examples of inadequate operation, design, and safety risks resulting from poor understating of the particle-particle interaction in wet conditions (see Table 1).In coal gasification and biomass thermal conversion, particle agglomeration induced by liquid (tar) binding is a major problem.This occurs in fluidized bed reactors when the bed temperature falls below ~300 • C, either due to poor mixing or heat losses.The same agglomeration phenomena may occur due to the partial melting of particles when operating above the melting point of the bed material.

Process Problems
Coal/biomass gasification Oil/tar leading to agglomeration and severe degradation in the low-temperature regions of fluidized bed reactors

Catalytic cracking
Surface catalyst melting at high temperatures leads to dead zones and de-fluidization.

Pneumatic conveying
Moisture leading to solid slugging, high-pressure drop, wear, and line blockage Fluidized bed coating/drying Liquid presence leads to undesired agglomeration and particle segregation Theoretical modeling and simulation of wet fluidization and fluidized bed reactors in particular are challenging due to the high complexity of the particle-particle and the wall-particle interactions.The few reported models in this case are mainly focused on discrete element modeling (DEM) and coupled CFD-DEM approaches.These have been used for the simulation of some relevant wet fluidization cases such as in spray granulation and fluidized bed coating [17,18].The DEM approach is based on Newton's second law of motion to track individual particles in the flow field and the inter-particle contact forces are calculated by a linear or non-linear spring and a dashpot model to account for the particle inelasticity.The gas phase is treated as a continuum phase (Eulerian) and therefore its flow field is commonly computed from the solution of volume-averaged Navier-Stokes equations, thus creating what is sometimes denoted by Eulerian-Lagrangian approach.While these studies demonstrate the predictive capabilities of the DEM and its complementary approach to practical experimental studies, it suffers from a serious drawback that makes it unattractive for the simulation of realistic processes where a large number of particles exist.The implementation of momentum transfer at the single-particle level, as used by DEM models, has a huge impact on the computational time.
The widely used two-fluid model or the Eulerian-Eulerian method is an alternative option to simulate solid-gas flow systems within a reasonable computational time at a realistic scale.However, the classic two-fluid based model, in its original formulations based on the KTGF, lacks the framework to incorporate the effect of the liquid layer, thus its predictive capabilities are limited to the dry flow case only.It is therefore widely recognized that it is important to identify new approaches for modeling fluidized flows at slightly wet conditions relevant to industrial applications.This research has also been given priority for future research focus at several recent meetings and conferences on granular flows and fluidization technology.In this study, the main objective of this study is to develop a continuum CFD model for the fluidization of slightly wet non-porous particles in non-agglomerating situations.First, the Electrical Capacitance Tomography (ECT) imaging technique is utilized to show the critical changes in the flow pattern and distinguishable flow regimes in a fluidized bed of glass beads wetted with small amounts of non-volatile liquid.This study then presents in detail the development and implementation of new constitutive equations to describe the particle-particle interaction in a slightly wet condition.To demonstrate the qualitative and quantitative validity of the model, the developed equations were formulated as a user-defined function (UDF) and implemented in a continuum CFD model for a bubbling fluidized bed using ANSYS Fluent CFD commercial code for multiphase flow.Comparisons between dry/wet fluidized bed experimental observation and the model predictions have been made.

Theory
In modeling gas-solid flow, knowledge of the solid shear stress is essential to describe the normal and tangential particle momentum transfer.In the case of slightly wet gas-solid flow, liquid bridges form between the contacting particles as demonstrated in Figure 1a.This results in changing the particle momentum due to inter-particle forces introduced by fluid shear resistance, while the collisional contacts dissipate energy in the liquid layer or bridges.Figure 1b demonstrates the basic hypothesis of the proposed flow regimes and the associated inter-particle forces.In the rapid flow regime, it is assumed that solid stress and dissipation are mainly dominated by random velocity fluctuations, and hence, the flow can be predicted using the original KTGF formulation without losing much accuracy.In the intermediate flow regime, the existence of a liquid layer surrounding the particles becomes more relevant; this is assumed to trigger a process in which the collision force normal to the particle is reduced.Energy is dissipated in the liquid layer and the tangential solid shear stress is substituted by a fluid shear resistance due to the formation of a liquid bridge.In the limiting regime of slow frictional flow, direct particle-particle enduring contacts take place due to the squeezing of the liquid bridges.In this case, it is assumed that the flow obeys the conventional dry frictional flow, which can be predicted by the critical state theory and onset of yielding as described by Coulomb law [19].
Energies 2024, 17, x FOR PEER REVIEW 4 of 20 The widely used two-fluid model or the Eulerian-Eulerian method is an alternative option to simulate solid-gas flow systems within a reasonable computational time at a realistic scale.However, the classic two-fluid based model, in its original formulations based on the KTGF, lacks the framework to incorporate the effect of the liquid layer, thus its predictive capabilities are limited to the dry flow case only.It is therefore widely recognized that it is important to identify new approaches for modeling fluidized flows at slightly wet conditions relevant to industrial applications.This research has also been given priority for future research focus at several recent meetings and conferences on granular flows and fluidization technology.In this study, the main objective of this study is to develop a continuum CFD model for the fluidization of slightly wet non-porous particles in non-agglomerating situations.First, the Electrical Capacitance Tomography (ECT) imaging technique is utilized to show the critical changes in the flow pattern and distinguishable flow regimes in a fluidized bed of glass beads wetted with small amounts of non-volatile liquid.This study then presents in detail the development and implementation of new constitutive equations to describe the particle-particle interaction in a slightly wet condition.To demonstrate the qualitative and quantitative validity of the model, the developed equations were formulated as a user-defined function (UDF) and implemented in a continuum CFD model for a bubbling fluidized bed using ANSYS Fluent CFD commercial code for multiphase flow.Comparisons between dry/wet fluidized bed experimental observation and the model predictions have been made.

Theory
In modeling gas-solid flow, knowledge of the solid shear stress is essential to describe the normal and tangential particle momentum transfer.In the case of slightly wet gas-solid flow, liquid bridges form between the contacting particles as demonstrated in Figure 1a.This results in changing the particle momentum due to inter-particle forces introduced by fluid shear resistance, while the collisional contacts dissipate energy in the liquid layer or bridges.Figure 1b demonstrates the basic hypothesis of the proposed flow regimes and the associated inter-particle forces.In the rapid flow regime, it is assumed that solid stress and dissipation are mainly dominated by random velocity fluctuations, and hence, the flow can be predicted using the original KTGF formulation without losing much accuracy.In the intermediate flow regime, the existence of a liquid layer surrounding the particles becomes more relevant; this is assumed to trigger a process in which the collision force normal to the particle is reduced.Energy is dissipated in the liquid layer and the tangential solid shear stress is substituted by a fluid shear resistance due to the formation of a liquid bridge.In the limiting regime of slow frictional flow, direct particleparticle enduring contacts take place due to the squeezing of the liquid bridges.In this case, it is assumed that the flow obeys the conventional dry frictional flow, which can be predicted by the critical state theory and onset of yielding as described by Coulomb law [19].

Solid Stress Tensor
The hypotheses described above have been used to derive the appropriate mathematical formulations required to take into consideration the reduced particle elasticity and increased inter-particle forces in particles of a slightly wet surface.Within the framework of continuum modeling, a strategy is proposed here to add a "liquid" normal stress term into the momentum equation of the particle phase.By analogy with the widely used kineticfrictional solid stress model, a modified solid shear stress tensor, taking into consideration the contribution of liquid-induced stresses, is defined by: where S is the strain rate, µ s,kin , µ s,col and µ s, f r represent the kinetic, collisional, and frictional contributions to the total solid-phase equivalent "viscosity", and these can be obtained from the existing classical KTGF and Coulomb model of friction.µ wet represents the viscosity arising from the particle-particle enduring contacts through the liquid bridge.Note that the model liquid used in this study is silicone oil.In this case, our calculation indicates that the capillary forces are insignificant compared to the viscous stresses, and hence ignored in the formulation of the solid stress model.It is also believed that the capillary force only becomes important in systems with high relative particle velocity or high-velocity fluctuation [20].
It is assumed that for a pair of wet contacting particles moving at a given velocity, .
u, the normal liquid force at the single-particle level can be defined by [21]: where µ L is the liquid dynamic viscosity, d s is the particle diameter and h is the inter-particle gap distance.The interparticle approach velocity can be estimated from the granular temperature, Θ s , as follows [22]: .
With the normal liquid bridge force at the particle level is in hand, further development is required to allow the implementation of this into continuum modeling framework, hence making the theory more appealing in the simulation of realistic wet solid suspension systems.In this case, the liquid bridge force must be expressed in terms of force per unit area (stress) within an ensemble.In doing so, the number of particles per unit area was first determined as follows: where α s is the solid concentration within an ensemble.Combining Equations ( 2)-( 4), a constrained normal solid stress is proposed as follows: where h is the inert-particle distance, h a is the characteristic measure of the particle surface asperity and h critical is the critical rupture distance of the liquid bridge.Note that Equation ( 5) ignores the effect of multi-bridges on particle, a simplification made to avoid using an adjustable parameter without strong physical justification.
Energies 2024, 17, 2656 6 of 20 By analog to Mohr-Coulomb law for a dry particle at enduring contact or yielding condition, an equivalent liquid shear viscosity is defined by: where ϕ ′ is the maximum angle of stability.Experiments on the stability of silicone oilwetted glass beads carried out by Nowak et al. [23] suggested a value ϕ ′ in the range of 20 • to 40 • .In this study, an average value of 30 • was used.The critical separation distance h critical , appearing in Equation ( 5), is defined here as the limit of the inter-particle gap between two contacting wet particles after which the liquid bridge rupture occurs.

Inter-Particle Gap and Liquid Bridge Rupture Distance
In randomly packed spheres, Woodcock [24] proposed a relation between the solid volume fraction and the inter-particle gap distance expressed in terms of the particle diameter as follows: In obtaining the critical rupture distance for particles at contact, h critical , a correlation equation derived from the original work of Lian et al. [25] has been used as follows: where γ is the liquid contact angle (in radians) and the term under the cube root represents the liquid-to-solid volume ratio.M L /M s is the overall ratio of liquid to solid mass in the bed (in the rest of this paper this ratio is denoted by δ). Figure 2a shows the relation between the separation distance and the solid volume fraction for the particle size considered in this study.It is shown that the separation distance decreases with increasing the solid concentration and correctly reduces towards zero as the solid volume fraction approaches the maximum packing condition.Figure 2b shows the critical rupture distance for the conditions considered here.The curve trend indicates that h critical is highly sensitive to the variation in liquid presence within the low range (δ < 0.03 × 10 −2 ), after which the variation becomes negligible.
Energies 2024, 17, x FOR PEER REVIEW 6 of 20 By analog to Mohr-Coulomb law for a dry particle at enduring contact or yielding condition, an equivalent liquid shear viscosity is defined by: where  is the maximum angle of stability.Experiments on the stability of silicone oilwetted glass beads carried out by Nowak et al. [23] suggested a value  in the range of 20° to 40°.In this study, an average value of 30° was used.The critical separation distance ℎ , appearing in Equation ( 5), is defined here as the limit of the inter-particle gap between two contacting wet particles after which the liquid bridge rupture occurs.

Inter-Particle Gap and Liquid Bridge Rupture Distance
In randomly packed spheres, Woodcock [24] proposed a relation between the solid volume fraction and the inter-particle gap distance expressed in terms of the particle diameter as follows: In obtaining the critical rupture distance for particles at contact, ℎ , a correlation equation derived from the original work of Lian et al. [25] has been used as follows: where  is the liquid contact angle (in radians) and the term under the cube root represents the liquid-to-solid volume ratio. / is the overall ratio of liquid to solid mass in the bed (in the rest of this paper this ratio is denoted by δ). Figure 2a shows the relation between the separation distance and the solid volume fraction for the particle size considered in this study.It is shown that the separation distance decreases with increasing the solid concentration and correctly reduces towards zero as the solid volume fraction approaches the maximum packing condition.Figure 2b shows the critical rupture distance for the conditions considered here.The curve trend indicates that ℎ is highly sensitive to the variation in liquid presence within the low range (δ < 0.03 × 10 −2 ), after which the variation becomes negligible.

Energy Dissipation
The KTGF assumes that energy dissipation arises from the inelastic collision of dry smooth particles.While this appears to be satisfactory in simulating various granular flow conditions, there is an argument that the dissipation term derived from the KTGF has shortcomings due to neglecting the dissipation caused by sliding friction [26].In the wet case, it is particularly important to take into consideration the increased energy dissipation due to non-ideal particle interactions (i.e., through a liquid bridge for example).In this proposed model, the energy dissipation induced by the liquid shear is implicitly considered through the change in the particle momentum equation, and the energy dissipation due to the normal collision is incorporated by using a modified particle restitution coefficient in the widely used energy dissipation equation derived from the KTGF as follows [27]: where e wet represents the restitution coefficient of wet particles.Davis et al. [28] measured the rebound velocities of small spherical particles falling on a smooth surface coated with a layer of viscous fluid at impact velocities within the range of 1-5 m/s.Accordingly, an expression for the wet particle restitution coefficient was derived from the Stokes analysis of the collision kinetic energy as follows: where e dry is the dry particle restitution coefficient, and St is the dimensionless Stokes number, which is a measure of the ratio of granular kinetic energy to the viscous dissipation and St * is the critical Stokes number, above which the liquid bridge rupture becomes important.To give allowance for variation of energy dissipation during liquid bridge deformation and complete loss of elasticity at enduring contacts, a constrained term for the restitution coefficient based on the Stokes analysis and the inter-particle gap was proposed as follows: Equation (11) allows for (i) a reasonable approximation of the energy dissipation in the liquid layer within the gap distance h a < h ≤ h critical (ii) complete dissipation when the particles collide with very low kinetic energy, and (iii) diminishing the effect of liquid layer dissipation while the particles are not subjected to liquid viscous stress within the gap distance of h > h critical .The Stokes number is defined as follows [29]: where m is the mass of the particle and µ L is the liquid dynamic viscosity.Thielmann et al.
[30] identified the critical Stokes number, St * , in terms of the inter-particle gap distance, h, and the dry particle restitution coefficient as follows: For smooth glass beads, h a falls within the range of 0.08-6 µm [6,30].In this study, an average value of h a = 2.8 µm has been used.

Experiments
Experiments were carried out in two different sized fluidization columns to collect quantitative and qualitative data on the hydrodynamics of slightly wet glass beads in a fluidized bed.The larger column was made of a cast acrylic column of 150 cm in height and 15 cm in diameter and equipped with a twin-plane ECT imaging system for recording the variations in solid concentration during fluidization.Air at ambient conditions was used to fluidize the bed material consisting of narrow-sized glass beads mixture of 350 µm mean diameter.Further details on this fluidization set-up, the ECT imaging system, and its application to fluidized bed hydrodynamic analysis and the image reconstruction procedure can be found in Makkawi and Wright [7,31] and Makkawi and Ocone [32].A few additional experiments were carried out in a smaller glass column of 100 cm height and 5 cm diameter.Both columns were fitted with a pressure probe located 1 cm above the gas inlet and connected to a pressure transducer and data recording system.To ensure reproducibility, all the experiments were repeated at least three times.
At the beginning of experiments, the fluidized bed was operated at fully turbulent conditions and a controlled amount of non-volatile silicone oil (Fluid 500, Dow Corning Ltd., Barry, UK) was added from the top using a syringe.Silicone oil is mainly used because it is non-volatile, hence allowing for mimicking wet fluidization in several practical conditions.Besides, silicone oil is of low permittivity, making it ideal for ECT measurement because it does not interfere with its signal.To ensure adequate wetting of particles, the bed was left to fluidize for at least 30 min.The air was then switched off and restarted again for fluidization recording at various gas velocities up to 0.8 m/s (this is 4.7 times the minimum fluidization velocity of the dry bed, i.e., U mf = 0.17 m/s).The pressure drop across the distributor in the fluidization columns was designed to give around 40% of the dry bed pressure drop, which falls within the "rule of thumb" for distributor design [33].In each experiment, the liquid loading was incrementally varied up to a maximum of δ = 0.15 × 10 −2 (mass of liquid to mass of dry bed).Additional details of the experiment operating conditions are given in Table 2.

Computational Model
The main model equation to describe the continuity, momentum and kinetic energy equations are described in this section as follows: Continuity equations: Energies 2024, 17, 2656 9 of 20 Momentum equations: The granular temperature was given by solving the pseudo kinetic energy equation given as follows [34]: The above model, with the constitutive equations given in Table 3 applies to the dry case.In simulating for the wet cases, the same set of the above main equations were solved but with the modified solid stress and energy dissipation, as given earlier in Section 3. The model was solved in three-dimensional coordinates using ANSYS Fluent CFD code (Ver 2016).Note that the effect of liquid presence on the kinetic energy (Equation ( 18)) is considered through the implementation of the modified solid stress and dissipation terms (implicitly in terms of modified restitution coefficient), as described in Section 3. (αs,max−αs) 3   Radial distribution function: Gas-solid drag coefficient: Diffusion coefficient of granular energy: 384(e ss +1)g 0,ss 1 + 6 5 α s g 0,ss e s,dry + 1 2 + 2α s 2 ρ s d s g 0,ss e s,dry Kinetic viscosity: 1 + 2 5 e s,dry + 1 3e s,dry − 1 α s g 0,ss Collisional viscosity: µ s,col = 4  5 α s ρ s d s g 0,ss e s,dry Frictional viscosity: Strain rate (i = gas or solid):

Meshing of the Simulation Domain and Solution Procedure
Figure 3 shows the geometry and meshing of the simulation domain using structured grids.The grid face size was set to minimum and maximum sizes of 3 and 7 mm, respectively, giving a total number of cells of 133,760.To confirm the reliability of meshing, grid-size sensitivity analysis was performed and this was found to produce an acceptable level of grid-independent solution.To ensure easy convergence of the partial differential equations (PDE) and optimize the computational time, the solution time step was determined using the Courant-Friedrichs-Lewy (CFL) condition, which is defined for three-dimensional PDE solution as follows: where C max is specified by the CFL condition to fall within the range of ~1-5 [35], ∆t is the solution time step and ∆x, ∆y and ∆z are the grid sizes in the x, y and z coordinates, respectively.u s,x , u s,y and u s,z are the particle velocity in the x, y and z coordinates, respectively.In this study, two different time steps of ∆t = 0.1 × 10 −3 and 0.3 × 10 −4 s were used for dry and wet conditions, respectively.The wet case was found to require a smaller time step to avoid the stiffness and instability in the solution of the solid momentum equation.The convection terms in the continuity, momentum and granular temperature equations were solved using first-order discretization schemes with a maximum number of iterations of 20 per time step and a convergence criterion of 10

Meshing of the Simulation Domain and Solution Procedure
Figure 3 shows the geometry and meshing of the simulation domain using structured grids.The grid face size was set to minimum and maximum sizes of 3 and 7 mm, respectively, giving a total number of cells of 133,760.To confirm the reliability of meshing, gridsize sensitivity analysis was performed and this was found to produce an acceptable level of grid-independent solution.To ensure easy convergence of the partial differential equations (PDE) and optimize the computational time, the solution time step was determined using the Courant-Friedrichs-Lewy (CFL) condition, which is defined for three-dimensional PDE solution as follows: where Cmax is specified by the CFL condition to fall within the range of ~1-5 [35], ∆t is the solution time step and ∆x, ∆y and ∆z are the grid sizes in the x, y and z coordinates,

Boundary Conditions and Simulation Parameter
The boundary condition for the gas phase at the wall was assumed to obey no-slip condition (zero velocity), while the solid velocity and granular temperature were assumed to satisfy the boundary condition proposed by Johnson and Jackson as follows [37]: (3Θ s,w ) 1/2 π φρ s α s g 0,ss δu s,w δn (20) where φ is the specularity coefficient and e wall is the particle-wall restitution coefficient.To simplify the problem and avoid further increase in the computational time without jeopardizing the solution accuracy, the particle-wall restitution coefficient was kept independent of the bed wetting condition.A summary of the parameters used in the model are given in Table 4.

Boundary Conditions and Simulation Parameter
The section described the experimental results obtained for the large fluidized bed using the ECT system.Figure 4 shows a stack of the ECT slice images taken during the fluidization of dry and wet glass beads.These images represent a "cross-sectional view" of the bed at a given level above the distributor.It is clear in Figure 4 that a small amount of liquid in a bubbling fluidized bed causes dramatic changes in the solid distribution and bubbles characteristics.These changes can be summarized as follows: i.
At a small liquid presence (δ < 0.08 × 10 −2 ) there is a decrease in the bubble size, as a result of bubbles splitting, and intensified wall-bubble movements.This may also change to gas channeling with most of the gas rising through a channel adjacent to the wall.In some way, this resembles the fluidization behavior of Geldart A dry particles [2].ii.
At an intermediate liquid presence (0.08 × 10 −2 < δ < 0.12 × 10 −2 ) slugging is observed (jointly rising bubbles occupying more than half of the column crosssection), typically causing the bed to rise and collapse as a piston.This is of great similarity to the fluidization behavior of highly cohesive powder such as Geldart C dry particles [2].iii.
At a considerably large liquid presence (δ > 0.12 × 10 −2 ) the bed becomes defluidized with the gas mainly escaping through one or more channels.
Quantitatively, the increase in bubbles splitting at a small liquid presence can be assessed by measuring the bubbles frequency.For this purpose, the time-series data of the solid fraction fluctuation recorded by the ECT system (representing the spatial average fluctuation at the height of 7.6 cm above the distributor) was analyzed using the Fast Fourier Transform (FFT) method.A comparison of the bubble frequency spectra for dry and wet beds is shown in Figure 5.It is demonstrated that the wet bed has a wider bubble frequency distribution and more than doubled dominant frequency due to the intensified bubble splitting, thus lending support to the observation noted in point (i) above.iii.At a considerably large liquid presence (δ > 0.12 × 10 −2 ) the bed becomes de-fluidized with the gas mainly escaping through one or more channels.Quantitatively, the increase in bubbles splitting at a small liquid presence can be assessed by measuring the bubbles frequency.For this purpose, the time-series data of the solid fraction fluctuation recorded by the ECT system (representing the spatial average fluctuation at the height of 7.6 cm above the distributor) was analyzed using the Fast Fourier Transform (FFT) method.A comparison of the bubble frequency spectra for dry and wet beds is shown in Figure 5.It is demonstrated that the wet bed has a wider bubble frequency distribution and more than doubled dominant frequency due to the intensified bubble splitting, thus lending support to the observation noted in point (i) above.iii.At a considerably large liquid presence (δ > 0.12 × 10 −2 ) the bed becomes de-fluidized with the gas mainly escaping through one or more channels.Quantitatively, the increase in bubbles splitting at a small liquid presence can be assessed by measuring the bubbles frequency.For this purpose, the time-series data of the solid fraction fluctuation recorded by the ECT system (representing the spatial average fluctuation at the height of 7.6 cm above the distributor) was analyzed using the Fast Fourier Transform (FFT) method.A comparison of the bubble frequency spectra for dry and wet beds is shown in Figure 5.It is demonstrated that the wet bed has a wider bubble frequency distribution and more than doubled dominant frequency due to the intensified bubble splitting, thus lending support to the observation noted in point (i) above.

Comparison of the Experiment Results with the Model Predictions
This section compares the model predictions with the experimental observation in the small fluidized bed in liquid presence.In Figure 6, it is evident that the model correctly predicts the four characteristic features of wet fluidization as discussed above.The bubbles splitting, gas channeling, and slugging associated with the liquid presence are demonstrated in Figure 6b-d.

Comparison of the Experiment Results with the Model Predictions
This section compares the model predictions with the experimental observation in the small fluidized bed in liquid presence.In Figure 6, it is evident that the model correctly predicts the four characteristic features of wet fluidization as discussed above.The bubbles splitting, gas channeling, and slugging associated with the liquid presence are demonstrated in Figure 6b-d.Quantitative validation is considered here by comparing the fluidization quality.Fluidization quality is characterized by a fluidization index (FI), commonly defined in terms of the bed pressure drop (∆P), the mass of the bed material (M bed ) and the bed cross-sectional area (A) as follows: Energies 2024, 17, 2656 14 of 20 Pressure drop measurement remains the simplest, yet the most effective method to detect changes in the flow pattern and identify the regime transition and minimum fluidization velocities in a fluidized bed.If the dry particles are fully suspended by the gas stream, then the FI is expected to be close to unity, since in this condition the pressure force ideally balances the bed weight force.Figure 7 shows the measured pressure drop at various liquid presence and gas velocities.It is interesting to note that, compared to the dry bed measurement, the pressure drop curve for the cases of the low liquid content of δ = 0.027 × 10 −2 and δ = 0.055 × 10 −2 is lower than that of the dry case and appears to show a peak at the point of incipient fluidization velocity of 0.2 m/s.In this range of wetting conditions, the pressure drop peak is believed to be associated with the breakdown of interparticle cohesion force and particle rearrangement before fluidization.The peak is then followed by reduced pressure due to the channeling phenomena.At the high range of liquid presence of δ > 0.1 × 10 −2 , the pressure drop is consistently higher than that of the dry case, mainly due to the commencement of slugging and de-fluidization.
Quantitative validation is considered here by comparing the fluidization quality.Fluidization quality is characterized by a fluidization index (), commonly defined in terms of the bed pressure drop (∆P), the mass of the bed material ( ) and the bed cross-sectional area () as follows: Pressure drop measurement remains the simplest, yet the most effective method to detect changes in the flow pattern and identify the regime transition and minimum fluidization velocities in a fluidized bed.If the dry particles are fully suspended by the gas stream, then the  is expected to be close to unity, since in this condition the pressure force ideally balances the bed weight force.Figure 7 shows the measured pressure drop at various liquid presence and gas velocities.It is interesting to note that, compared to the dry bed measurement, the pressure drop curve for the cases of the low liquid content of δ = 0.027 × 10 −2 and δ = 0.055 × 10 −2 is lower than that of the dry case and appears to show a peak at the point of incipient fluidization velocity of 0.2 m/s.In this range of wetting conditions, the pressure drop peak is believed to be associated with the breakdown of interparticle cohesion force and particle rearrangement before fluidization.The peak is then followed by reduced pressure due to the channeling phenomena.At the high range of liquid presence of δ > 0.1 × 10 −2 , the pressure drop is consistently higher than that of the dry case, mainly due to the commencement of slugging and de-fluidization.The comparison of measured and predicted  is shown in Figure 8.As expected, the  value at dry conditions is close to unity.The overall model-based predictions and experimental curve trends are in satisfactory agreement within the low range of liquid loading.The observed decrease in  is mainly due to the increase in bubbles splitting and channelling at this limit of wetting condition, as demonstrated earlier in Figure 6.At a higher liquid loading of δ > 0.1 × 10 −2 , the experimental data shows the value of  to increase above unity due to severe slugging (piston-like behavior) and increased wall The comparison of measured and predicted FI is shown in Figure 8.As expected, the FI value at dry conditions is close to unity.The overall model-based predictions and experimental curve trends are in satisfactory agreement within the low range of liquid loading.The observed decrease in FI is mainly due to the increase in bubbles splitting and channelling at this limit of wetting condition, as demonstrated earlier in Figure 6.At a higher liquid loading of δ > 0.1 × 10 −2 , the experimental data shows the value of FI to increase above unity due to severe slugging (piston-like behavior) and increased wall shearing and friction, eventually leading to de-fluidization.This is to a great extent similar to the behavior of a highly cohesive powder, where the FI value is usually greater than 1.4 [38].The model fails to provide a solution at a high liquid loading (δ > 0.1 × 10 −2 ) because at this condition the bed is virtually static (de-fluidized) and hence, the solution diverges.
shearing and friction, eventually leading to de-fluidization.This is to a great extent similar to the behavior of a highly cohesive powder, where the  value is usually greater than 1.4 [38].The model fails to provide a solution at a high liquid loading (δ > 0.1 × 10 −2 ) because at this condition the bed is virtually static (de-fluidized) and hence, the solution diverges.

Solid Shear Stress, Energy Dissipation, and Granular Temperature
The numerical values of the dry flow kinetic, collisional, and frictional shear stress at different flow regimes (rapid, intermediate, and dense) have been previously reported by Makkawi and Ocone [39] and Makkawi et al. [40] using a one-dimensional two-fluid model based on the kinetic theory of granular flow.It was shown that at high solid fraction, the shear stress approaches the Coulomb yield condition, while in dilute flow, the particle stresses are dominated by the kinetic contribution.In between these boundaries, there is a mixed contribution of collisional and kinetic contributions.In a slightly wet particle flow, the model proposed here assumes that the liquid-induced bridge shear forces mainly dominate the flow at the intermediate regime, which has its boundaries identified in terms of the solid volume fraction through the interparticle gap distance.
Figure 9 shows the predicted solid shear stresses in the fluidized bed at the liquid loading of δ > 0.1 × 10 −2 .Due to the anisotropic nature of the stresses, the reported values are for the root sum square given by  +  +  ).The values of the kinetic, collisional, and frictional contribution to the overall shear stress are numerically close to the range reported in the literature [41].The liquid-induced shear stress within the intermediate flow regime dominates the overall shear stress and approaches the frictional stress at close to maximum packing.It should be noted that the scattering of the various solid stress terms appearing in this figure, i.e., different stress values at the same volume fraction, is owing to the variations of the shear rate and granular temperature at the same solid concentration.The numerical values of the dry flow kinetic, collisional, and frictional shear stress at different flow regimes (rapid, intermediate, and dense) have been previously reported by Makkawi and Ocone [39] and Makkawi et al. [40] using a one-dimensional two-fluid model based on the kinetic theory of granular flow.It was shown that at high solid fraction, the shear stress approaches the Coulomb yield condition, while in dilute flow, the particle stresses are dominated by the kinetic contribution.In between these boundaries, there is a mixed contribution of collisional and kinetic contributions.In a slightly wet particle flow, the model proposed here assumes that the liquid-induced bridge shear forces mainly dominate the flow at the intermediate regime, which has its boundaries identified in terms of the solid volume fraction through the interparticle gap distance.
Figure 9 shows the predicted solid shear stresses in the fluidized bed at the liquid loading of δ > 0.1 × 10 −2 .Due to the anisotropic nature of the stresses, the reported values are for the root sum square given by τ 2 xy + τ 2 xz + τ 2 yz ).The values of the kinetic, collisional, and frictional contribution to the overall shear stress are numerically close to the range reported in the literature [41].The liquid-induced shear stress within the intermediate flow regime dominates the overall shear stress and approaches the frictional stress at close to maximum packing.It should be noted that the scattering of the various solid stress terms appearing in this figure, i.e., different stress values at the same volume fraction, is owing to the variations of the shear rate and granular temperature at the same solid concentration.Figure 10 shows the dependence of the energy dissipation and granular temperature on the solid concentration.In Figure 10a, the energy dissipation for the wet case is much higher than for the dry case.This is due to the reduced particle elasticity and dissipation in the liquid bridge as the restitution coefficient drops to low values.In Figure 10b, the granular temperature appears to generally drop with increasing the solid concentration.This is expected since as the particle approaches dense packing conditions, little room is available for relative particle motion and collision.Similar trends and numerical values for the dry case have been reported by Gidaspow et al. [41] and Makkawi and Ocone [39].In the wet condition, the granular temperature is expected to be lower due to the lower particle velocity fluctuation, however, this is not well depicted here.It is possible that the energy generation associated with the additional liquid viscous stress, as implemented in the energy equation (Equation ( 18)), contributes to increasing the granular temperature, which is then counterbalanced by the increased energy dissipation as evident in Figure 10a.In such a case, the granular temperature may well remain within the range of the dry case due to the competing dissipation and generation terms.We believe that the relation between the energy dissipation and granular temperature, particularly for the wet case, is a complex one and requires further investigation.Figure 10 shows the dependence of the energy dissipation and granular temperature on the solid concentration.In Figure 10a, the energy dissipation for the wet case is much higher than for the dry case.This is due to the reduced particle elasticity and dissipation in the liquid bridge as the restitution coefficient drops to low values.In Figure 10b, the granular temperature appears to generally drop with increasing the solid concentration.This is expected since as the particle approaches dense packing conditions, little room is available for relative particle motion and collision.Similar trends and numerical values for the dry case have been reported by Gidaspow et al. [41] and Makkawi and Ocone [39].In the wet condition, the granular temperature is expected to be lower due to the lower particle velocity fluctuation, however, this is not well depicted here.It is possible that the energy generation associated with the additional liquid viscous stress, as implemented in the energy equation (Equation ( 18)), contributes to increasing the granular temperature, which is then counterbalanced by the increased energy dissipation as evident in Figure 10a.In such a case, the granular temperature may well remain within the range of the dry case due to the competing dissipation and generation terms.We believe that the relation between the energy dissipation and granular temperature, particularly for the wet case, is a complex one and requires further investigation.

Conclusions
The proposed model allowed, for the first time, continuum modeling of slightly wet solid-gas fluidized bed, thus extending the existing two-fluid modeling beyond its traditional boundaries.The proposed model combines theories of liquid bridge forces with the kinetic theory of granular flow (KTGF) formulations.The significant impact of the liquid presence is depicted by the increased energy dissipation similar to that observed in highly cohesive powder.The liquid-induced shear stress is found to dominate the overall shear

Conclusions
The proposed model allowed, for the first time, continuum modeling of slightly wet solid-gas fluidized bed, thus extending the existing two-fluid modeling beyond its traditional boundaries.The proposed model combines theories of liquid bridge forces with the kinetic theory of granular flow (KTGF) formulations.The significant impact of the liquid presence is depicted by the increased energy dissipation similar to that observed in highly cohesive powder.The liquid-induced shear stress is found to dominate the overall shear stress and approaches the frictional stress at close to maximum packing.The experimental measurement produced by the electrical capacitance tomography (ECT) showed distinct hydrodynamic features characterized by bubbles splitting, gas channeling, slugging, and de-fluidization as the liquid present in the bed increases.Such experimental data is scarce in the open literature.Further experimental work at wider operating conditions (e.g., particle sizes, fluidization velocity, liquid properties) is highly recommended to reveal more details, provide important data for the validation of theoretical models, and help in the safe design and operation of modern fluidized bed systems for power generation and beyond.Finally, it is hoped that the proposed model will open the door for further development to address the increasing demand for new generations of fast computational models that can reasonably predict the flow in such a complex flow condition.
Author Contributions: Y.M. developed the research conceptualization, derived the model equations, developed the methodology of the solution, conducted the experimental investigation, led the project supervision and administration, and contributed to the data interpretation and writing/editing of the manuscript.X.Y. carried out the computational work using a software program (Fluent), implemented the constitutive equations on the software, and contributed to the data interpretation and writing of the draft manuscript.R.O. contributed to the research proposal, funding acquisition, project administration, and project planning.S.G. contributed to the research proposal, supervision, and writing of the manuscript.All authors have read and agreed to the published version of the manuscript.
Funding: This research was funded by The Leverhulme Trust (Grant: RPG-410).
Data Availability Statement: Data are available upon request.

Conflicts of Interest:
The authors declare no conflicts of interest.

C
Courant number (-) C D drag coefficient (-) d s particle diameter of solid phase (m) e dry dry particle-particle restitution coefficient (-) e ss particle-particle restitution coefficient (-) e s i ,w particle-wall restitution coefficient (-) e wet wet particle-particle restitution coefficient (-)

Figure 1 .
Figure 1.Wet particle interaction: (a) microscope image of wet particles (glass beads wetted with silicone oil) at static conditions demonstrating the formation of liquid bridges at the contact point; (b) schematics of the proposed flow regime and the corresponding particle-particle contacts.

Figure 2 .
Figure 2. (a) Inter-particle distance and (b) the critical separation distance for a particle diameter of 350 µm.

Figure 2 .
Figure 2. (a) Inter-particle distance and (b) the critical separation distance for a particle diameter of 350 µm.

Figure 3 .
Figure 3. Computation domain and meshing of the fluidization column.

Figure 3 .
Figure 3. Computation domain and meshing of the fluidization column.

Figure 4 .
Figure 4. Electric capacitance tomography (ECT) images of 7.6 cm above the fluidized bed gas distributor.Each column is produced by stacking a series of images recorded at the rate of 100 frames per second and representing the spatial average solid concentration over 3.8 cm height.The data was produced in a 15 cm diameter column with the bed material consisting of 3.5 kg glass beads fluidized at the gas velocity of 0.8 m/s.

Figure 5 .
Figure 5. Fast Fourier transform (FFT) analysis of solid fraction fluctuation obtained by ECT measurement in a bubbling fluidized bed (a) dry bed (b) wet bed at δ = 0.055 × 10 −2 .The data were produced at 0.8 m/s gas velocity.

Figure 4 .
Figure 4. Electric capacitance tomography (ECT) images of 7.6 cm above the fluidized bed gas distributor.Each column is produced by stacking a series of images recorded at the rate of 100 frames per second and representing the spatial average solid concentration over 3.8 cm height.The data was produced in a 15 cm diameter column with the bed material consisting of 3.5 kg glass beads fluidized at the gas velocity of 0.8 m/s.

Figure 4 .
Figure 4. Electric capacitance tomography (ECT) images of 7.6 cm above the fluidized bed gas distributor.Each column is produced by stacking a series of images recorded at the rate of 100 frames per second and representing the spatial average solid concentration over 3.8 cm height.The data was produced in a 15 cm diameter column with the bed material consisting of 3.5 kg glass beads fluidized at the gas velocity of 0.8 m/s.

Figure 5 .
Figure 5. Fast Fourier transform (FFT) analysis of solid fraction fluctuation obtained by ECT measurement in a bubbling fluidized bed (a) dry bed (b) wet bed at δ = 0.055 × 10 −2 .The data were produced at 0.8 m/s gas velocity.

Figure 5 .
Figure 5. Fast Fourier transform (FFT) analysis of solid fraction fluctuation obtained by ECT measurement in a bubbling fluidized bed (a) dry bed (b) wet bed at δ = 0.055 × 10 −2 .The data were produced at 0.8 m/s gas velocity.

Figure 6 .
Figure 6.Comparison of experimental and predicted fluidized bed material distribution at various wetting conditions.The experiments were carried out in a small fluidization column of 5 cm diameter at a fluidization velocity of 0.8 m/s.

Figure 6 .
Figure 6.Comparison of experimental and predicted fluidized bed material distribution at various wetting conditions.The experiments were carried out in a small fluidization column of 5 cm diameter at a fluidization velocity of 0.8 m/s.

Figure 7 .
Figure 7. Experimental fluidized bed pressure drop at various gas velocities and wetting conditions.The data was produced in a 15 cm diameter column with the bed material consisting of 3.5 kg glass beads.Each data point represents the average of three measurements (maximum deviation of ∓10%).

Figure 7 .
Figure 7. Experimental fluidized bed pressure drop at various gas velocities and wetting conditions.The data was produced in a 15 cm diameter column with the bed material consisting of 3.5 kg glass beads.Each data point represents the average of three measurements (maximum deviation of ∓10%).

Figure 8 .
Figure 8.Comparison of the predicted and experimentally determined fluidization index (FI).The experimental data was produced in a 15 cm diameter column with the bed material consisting of 3.5 kg glass beads fluidized at the gas velocity of 0.8 m/s.Each data point represents the average of three measurements (maximum deviation of ∓10%).

Figure 8 .
Figure 8.Comparison of the predicted and experimentally determined fluidization index (FI).The experimental data was produced in a 15 cm diameter column with the bed material consisting of 3.5 kg glass beads fluidized at the gas velocity of 0.8 m/s.Each data point represents the average of three measurements (maximum deviation of ∓10%).

Figure 9 .
Figure 9. Predicted solid shear stress in a slightly wet fluidized bed of 15 cm diameter at the gas velocity of 0.8 m/s and liquid presence of δ = 0.1 × 10 −2 wt% liquid.

Figure 9 .
Figure 9. Predicted solid shear stress in a slightly wet fluidized bed of 15 cm diameter at the gas velocity of 0.8 m/s and liquid presence of δ = 0.1 × 10 −2 wt% liquid.

Figure 10 .
Figure 10.Predicted (a) energy dissipation rate and (b) granular temperature as a function of the solid concentration in a dry and a slightly wet fluidized bed of 15 cm diameter at the gas velocity of 0.8 m/s.

Figure 10 .
Figure 10.Predicted (a) energy dissipation rate and (b) granular temperature as a function of the solid concentration in a dry and a slightly wet fluidized bed of 15 cm diameter at the gas velocity of 0.8 m/s.

Table 1 .
Examples of problems and risks associated with wet particle flow and suspension.

Table 3 .
Constitutive relations for the gas-solid flow.Solids pressure: P s = P s,kin + P s,col + P s, f r where P s,kin = α s ρ s Θ s P s,col = 2ρ s 1 + e s,dry α 2 s g 0,ss Θ s P s, f r = 0.05 (αs−α s,min ) 2

Table 4 .
Summary of the parameters used in the model.