Improvement of voidage prediction in liquid-solid fluidized beds by inclusion of the Froude number in effective drag relations

A novel effective drag relation for liquid-solid ﬂuidisation is proposed, suitable for application in full-scale installations. This is achieved by presenting new insights related to the inﬂuence of the temporal-spatial heterogeneity on the effective hydrodynamic drag for large ﬂuidised systems. While heterogeneous ﬂow behaviour can be predicted increasingly accurately in CFD simulations that explicitly model the heterogeneous solids distribution, for the operation of many large-scale applications it is infeasible to perform such computationally intensive simulations. Therefore, there is a clear need for full-scale drag relations that effectively take into account the heterogeneous behaviour and irregular spatial particle distributions. Our new drag relation is based on a large set of experiments, which shows that the degree of overall expansion is not only dependent on the ratio of laminar-turbulent ﬂow, but also on the amount of homogenous versus heterogeneous ﬂow, which is not included in current full-scale drag relations. To include the effect of heterogeneity, the standard drag relation, based on the Reynolds number, is extended with a speciﬁc type of Froude number. Because fully turbulent ﬂow regimes are rare in applications of liquid-solid ﬂuidisation, our focus is not on the turbulent ﬂow regime but instead on laminar and transi- tional ﬂow regimes. In these regimes, three types of models are investigated. The ﬁrst type is based on a theoretical similarity with terminal settling, the second is based on the semi-empirical Carman-Kozeny model, and the third is based on empirical equations using symbolic regression techniques. For all three types of models, coeﬃcients are calibrated on experimental data with monodisperse and almost spherical glass beads. The models are validated with a series of calcium carbonate grains applied in drinking water treatment processes as well as data obtained from the literature. Using these models, we show that the voidage prediction average relative error decreases from approximately 5% (according to the best literature equations which use Reynolds number only) to 1–2% (using both Reynolds and Froude number). This implies that our new models are more suitable for operational control in full-scale ﬂuidised bed applications, such as pellet softening in drinking water treatment processes.


Introduction
Liquid-solid fluidisation is frequently used in drinking water treatment processes, for instance in seeded crystallisation softening processes ( Crittenden et al., 2012 ). For optimal process conditions, i.e. fast calcium carbonate crystallisation ( van Schagen, 2009 ), a large surface area in the fluidised bed is required. Therefore, in the field of water treatment, and specifically in pellet softening liquid-solid fluidisation reactors, accurate voidage prediction models are crucial . In a fluidised bed, the voidage (or bed voidage), particle size, and physical properties of the fluid and particles are inter-related. In the literature, numerous multiphase models are given to predict the overall voidage in fluidised bed reactors, mainly aimed at gas-solid systems ( Di Felice, 1995 ;Gibilaro, 2001 ;Yang, 2003 ;Crowe and Group, 2006 ;Rhodes, 2008 ;Seville and Yu, 2016 ;Yates and Lettieri, 2016 ;Johnson, 2016 ). Specifically for liquid-solid systems, the overall voidage can be predicted using classical models, such as the Richardson-Zaki approach based on terminal settling velocity, or an improved version also using the incipient fluidisation . The other frequently applied method is based on the idea of a flow through an assumed collection of channels in a bed of particles ( Kozeny, 1927 ;Carman, 1937 ). The well-known Carman-Kozeny equation is an important drag relation for the determination of permeability in porous media, such as filters, as well as for the estimation of the voidage in a fluidised bed in water treatment processes. The simplicity of this model is a consequence of an evident direct relationship between the particle Reynolds number and the drag exerted by the fluid on the particles ( Bird et al., 2007 ) where viscous and inertial forces are balanced. This drag model, however, does not consider homogeneous or heterogeneous flow patterns. Homogeneous fluidisation for uniform particles was observed by Wilhelm and Kwauk (1948) and Richardson and Zaki (1954) . Fluidisation quality was characterised by ( Geldart, 1973 ) who made a distinction between non-bubbling (particulate, homogeneous) and bubbling (aggregative heterogeneous) fluidisation, however this was for gas-liquid fluidisation. According to Couderc (1985), Davidson et al., 1985 , two modes of fluidisation occur in liquid-solid fluidised systems, i.e. particulate and aggregative fluidisation ( Davidson et al., 1985 ). Based on experimental studies ( Didwania and Homsy, 1981 ;Ham et al., 1990 ;Li and Kwauk, 1994 ) at least five flow regimes can be identified: stable uniform fluidisation, particulate regime, wavy flow, wavy flow with traverse structure, fine scale turbulent flow, and bubbling regime as the liquid velocity is increased. Di Felice (1995) reported inhomogeneities in liquid-solid fluidisation systems where the degree of heterogeneous behaviour depended on particular system characteristics such as morphological particle properties, particle size distribution, and fluid-to-solid density ratio. Small particles with a density closer to that of the fluidising medium are more easily fluidised compared to large and heavy particles, since the gravitational pull is larger for the latter particles. The interparticle forces on small particles are relatively more important than the same forces acting on large particles, causing small particles to exhibit a certain velocity range of homogeneous expansion ( Beetstra, 2005 ). For larger particles, throughout the bed, large inhomogeneities may occur: these include clustering of particles and voids of water as well as velocity fluctuations due to bubbling and spouting effects.
The importance of the fluid-particle interaction problem is considerable ( Wu and Yang, 2019 ). The modelling of full-scale fluidised bed reactors is challenging because of their complex flow behaviour and numerous particle interactions ( Cornelissen et al., 2007 ). Using Direct Numerical Simulation (DNS]), the whole range of spatial and temporal scales can be resolved in the computational mesh, which therefore necessarily contains only a few particles (N~100 to 10 0 0). The advantage of DNS is that no explicit drag relations need to be imposed. Rather, the drag is resolved for each particle, taking into account the effects of particle-particle collisions and other forces acting on the particles ( Al-Arkawazi et al., 2017 ). The computational costs however are very high. The costs can be lowered by using a coupled Computational Fluid Dynamics -Discrete Element Method (CFD-DEM) model ( Ghatage et al., 2014 ), in which the mesh is actually larger than the particles, and therefore flow around the particles is no longer resolved. Although this allows for more particles to be simulated (N~10,0 0 0 to 10 0 0,0 0 0), a drawback is that the interactions between the fluid phase and solids phase must be modelled explicitly by a drag relation, for instance the Ergun model ( Liu et al., 2015 ). In CFD-DEM, it is important to not choose the mesh size too large either, to be able to assume a more-or-less homogeneous particle distribution at the scale of a single CFD cell. Regarding full-scale industrial fluidised bed reactors, like water pellet-softening (N~10,0 0 0,0 0 0,0 0 0), CFD simulations could perhaps be achievable, but only by approximating the particle interactions and particle-fluid drag even more, such as in Filtered Two-Fluid Models ( Cloete et al., 2018 ). Such CFD models are computationally very expensive, making them less suitable for process optimisation and plant wide control. To cope with constantly changing operational conditions in full-scale installations, more straightforward models are needed for optimal and robust process control. In particular, there is a need for drag re-lations that can predict the pressure drop and overall voidage in fluidised beds, effectively taking into account the local and global multiphase flow phenomena occurring in full-scale installations.
The aim of the current work is to determine an improved fullscale drag relation considering the abovementioned homogeneous and heterogeneous flow regimes. The majority of drag relations given in the literature is based on the Reynolds number only and assumes static and homogeneous particle arrangements. Because these drag relations ignore temporal-spatial fluctuations, they are not suitable to accurately cope with the transition from homogeneous to heterogenous flow regimes. Here, we hypothesise that the popular Carman-Kozeny model can be improved through the introduction of the Froude number, in addition to the Reynolds number, to incorporate the heterogeneous aspects. In fluid mechanics ( Wilkes, 2019 ), the dimensionless Froude number, is used to indicate the influence of gravity on fluid motion ( Di Felice, 1995 ). Wilhelm and Kwauk (1948) observed that the Froude number was a reliable parameter for discriminating between two extreme situations, with the fluidisation behaviour ranging from extremely smooth to violently bubbling. We will validate our new drag models by comparing their predictions with a large set of experimental results from our laboratories, as well as experimental data from the available literature.

Theory
Section 2.1 below elucidates the fundamental principles of drag relations, in particular Carman-Kozeny. The extension of the drag relation with the Froude number to cope with the heterogeneous flow aspects is introduced in Section 2.2. The nomenclature is given at the end of the manuscript.

Laminar flow regime: Blake-Kozeny equation
At low Reynolds numbers, the relation between the fluid flow velocity through a dense porous medium and the pressure drop over this system is described in general by Darcy's law. Accordingly, the head loss, or hydraulic gradient, has the following form for laminar flow in a packed bed of particles: The hydraulic gradient ( Ergun, 1952 ) is often given in terms of a drag coefficient or f T The dimensionless drag coefficient in Eq.
(2) is often given as a function of the modified particle Reynolds number Re ɛ : where the modified particle Reynolds number Re ɛ is defined as:

Laminar flow regime: Kozeny equation
The corresponding Kozeny drag coefficient f T is given as Kozeny (1927) :

Turbulent flow regime: Burke-Plummer equation
For the complete turbulent flow regime given by Burke and Plummer (1928) , the corresponding Burke-Plummer drag coefficient f T ( Bird et al., 2007 ) states: 2.1.4. Transitional flow regime: Ergun and Carman-Kozeny Ergun (1952) combined the Carman and Burke-Plummer equations and added these together, producing a mathematically blended model ( Eq. (7) ) to predict laminar, transitional, and turbulent flow, which satisfies the linear and nonlinear terms in the Reynolds number. The Ergun drag coefficient is based on experimental data between 2 < Re ɛ < 40 0 0 and is often used in CFD modelling ( Beetstra, 2005 ;Kolev, 2012 ;Erdim et al., 2015 ;Tavassoli et al., 2015 ): However, we will show below that there is a considerable discrepancy between the Ergun model and our experimental data for Re ɛ > 500.
For the transitional flow regime, the Carman-Kozeny drag coefficient ( Carman, 1937 ) is given by: The Carman-Kozeny equation is, de facto , the most commonly used equation and applied in various fields such as ground water flow, water treatment processes, and a variety of chemical engineering applications Xu and Yu, 2008 ). According to Sobieski and Zhang (2014) , both Ergun ( (7) and  equations are sensitive to voidages and more pronounced at lower voidage. Since in water treatment the operational field lies in the vicinity of incipient fluidisation and, in addition, turbulent flow regimes are exceptional, there is a preference for using the Carman-Kozeny drag relation. Additionally, van Dijk and Wilms (1991) proposed an empirical simplified drag relation based on the Kozeny Eq. (5) valid for the transition region.

Traditional drag versus Reynolds number
In general, drag relations, where the drag is given as a function of the Reynolds number, are commonly plotted on log-log scales to comprise a wide range of Reynolds numbers. The most common drag relation is pipe flow friction given in Moody's diagram ( Moody, 1944 ). Here, the laminar flow regime declines fast for increasing Reynolds numbers, as the emphasis lies on turbulent flow regimes, represented by horizontal lines.
For porous media, drag relations and plots in the literature show a similar pattern. The laminar drag declines with increasing Reynolds numbers following Kozeny ( Eq. (5) ), after which the drag deflects in the transitional flow regime towards the constant drag predicted by Burke-Plummer ( Eq. (6) ), which is valid in the turbulent flow regime ( Forchheimer, 1930 ). As turbulent flow regimes are rare in applications of liquid-solid fluidisation, the emphasis should not lie on the turbulent flow regime. To place greater emphasis on the drag in the laminar and transitional regions, it is better to focus on deviations from the laminar regime by multiplying f T (drag coefficient in the turbulent ( ~v s 2 ) representation) with the modified particle Reynolds number Re ɛ ( Eq. (4) ) to arrive at f L (drag coefficient in the laminar ( ~v s ) representation) ( Erdim et al., 2015 ): In the literature, deviations between measurements and drag models are commonly present, but often artificially hidden due to the use of log-log scales over many orders of magnitude. In this work, by reporting f L ( Eq. (11) ) values on a linear scale, we will make these deviations between measured and modelled drag more clearly visible in the transitional regime.

Drag relations for fixed and fluidised beds
Eq.
(2) can be rewritten to obtain an explicit expression for f T based on measurements of the hydraulic gradient: Another equation for the drag coefficient, denoted by f L , was given by Erdim et al. (2015) : In a steady state of homogeneous fluidisation of particulate solids, the frictional pressure drop equals the weight of the bed material, reduced by the buoyancy forces, per unit of bed surface. This is expressed by Eq. (14) ( Yang, 2003 ): Based on Eq. (14) , the force balance between the frictional drag and the weight of the particles in the bed yields an alternative expression for the dimensionless drag coefficient f T or f L : Eqs. (12) and ( 13 ) can be used to determine the dimensionless drag coefficient based on experimental data for both fixed and fluidised bed states, where the differential pressure and voidage must be known. Eqs. (15) and ( 16 ) are dependent on the voidage, albeit independent of the differential pressure, and they are valid only for the fluidised bed state. For accurate drag determination, Eq. (16) might be preferable since there is no dependency on sensitive differential pressure measurements. The determination of bed voidage in liquid-solid fluidisation is relatively accurate due to the straightforward measurement of total mass of particles, bed height, and particle density. Wilhelm and Kwauk (1948) proposed a dimensionless number, known as the Froude number, which is a good parameter for discriminating between the two extreme situations: particulate and bubbling behaviour, to explain the quality of fluidisation:

Hydraulic models based on the Reynolds and Froude numbers 2.2.1. Froude numbers
The Froude number can be viewed as the ratio of inertial to gravity forces. A transition occurs from particulate or smooth homogeneous fluidisation to heterogeneous or aggregative (bubbling) fluidisation at Fr ≈ 1. According to Gupta and Sathiyamoorthy (1999) , this transition occurs at Fr ≈ 0.13. In general, the transition from a homogeneous to a heterogeneous state is gradual. Many standard works define the Froude number as Eq. (17) , ( Di Felice, 1995 ;Gupta and Sathiyamoorthy, 1999 ;Bird et al., 2007 ;Yates and Lettieri, 2016 ;Rapp, 2017 ;Wilkes, 2019 ). Other standard works ( Yang, 2003 ;Rhodes, 2008 ) use the square root of the expression given in Eq. (17) : In case the particle and fluid densities are taken into account ( Grace, 1986 ;Crowe and Group, 2006 ), the densimetric or particle Froude number is given by: The densimetric Froude number, Eq. (19) , also contains information about the ratio of the densities of the particle and the fluid. Because the gravitational force on a particle in a fluid is always counteracted by a buoyancy force (a neutrally buoyant particle does not sediment or fluidise), it appears to be the more relevant Froude number.

Proposed model extensions
There are several reasons why we propose an extension of the existing classical drag models with the particle Froude number: -While the Reynolds number deals with the relationship between viscous and inertial forces, the particle Froude number deals with the relationship between gravity and inertial forces. Notably for larger particles and/or with a relatively large solid-to-fluid density, the gravitational forces are dominant. In other words: the Reynolds number quantifies the laminar-toturbulent properties, whilst the Froude number quantifies the homogeneous-to-heterogeneous properties of the system. -When aggregative fluidisation occurs, voids of fluid provide pathways with less resistance to the fluid, resulting in a lower drag compared to the drag assumed for homogeneous fluidisation. This means that, when we look at the form of equations (7) or (8) , the effective Reynolds number should in fact be slightly higher. This might be accomplished through adding an explicit dependency on the Froude number. -Visual observations of assumed homogeneous fluidisation experiments showed significant voids and trains of particles, which means that the drag cannot be mathematically described merely by the Reynolds number based on viscous and inertial forces. Through inclusion of the Froude number, information regarding such transient fluidisation events is appended. Videos with liquid-solid fluidisation experiments are shared in the Supplementary Material ( Kramer, 2019 ). -In DNS modelling ( Beetstra et al., 2007 ), drag relations are often proposed for static arrays of particles, whereas in many practical applications, such as pellet softening, the particles are moving in space, leading to heterogeneities as mentioned above.

Model synthesis incorporating the Reynolds and Froude numbers
To take into account both laminar and turbulent characteristics as well as homogeneous and heterogeneous phenomena, we propose three types of models for liquid-solid fluidised beds, based on the Reynolds number ( Eq. (4) ) along with the Froude number ( Eq. (19) ): The first, of our three models, is based on a theoretical similarity with terminal settling, the second is based on the semiempirical Carman-Kozeny model, and the third is based on empirical equations using symbolic regression techniques.
Our first model is based on theoretical expressions for the drag on a sphere falling through a quiescent fluid at small but finite Reynolds numbers. According to the principles introduced by Oseen (1927) and Proudman and Pearson (1957) , the singular solution of the basic Stokes equation ( Batchelor, 2012 ) grows likelog(Re) instead of decaying like 1/Re towards the transitional region. For increasing Reynolds numbers, the inertial forces become more important, and the drag coefficient approaches Newton's law i.e. a constant drag coefficient. The newly proposed dimensionless drag coefficient is a combination of the Reynolds number Re ɛ , and the Froude number With the application of our newly proposed dimensionless number RF , the form of Eq. (8) is therefore adjusted to: Previously, the numerical coefficients c 1 , c 2 , and c 3 in Eq. (23) were obtained by fitting them to experimental results ( Gupta and Sathiyamoorthy, 1999 ). In this work, however, the coefficients have been fitted through non-linear curve fitting. Note that for small Froude numbers (homogeneous fluidisation), RF should approach the modified Reynolds number Re ɛ . Still, this leaves many possibilities with regard to the exact relation in Eq. (22) . We have explored the following models: RIO 1 model (Reynolds-Improved-Outlook model) The first model is inspired by the shape of the equation proposed by Schiller and Naumann (1933) : The third model was obtained using symbolic regression techniques as applied in genetic programming. Genetic programming is a random-based technique ( Koza, 1992 ) for automatically learning computer programmes based on artificial evolution. It has been successfully used in many applications ( Edwards, 2006 ;Barati et al., 2014 ). The advantage of genetic programming is that there is no need to define the structure of a model a priori : the technique randomly generates a population of several mathematical operators. Symbolic regression is the process of determining the symbolic function, which describes a data set, thus effectively developing an empirical model ( Awange and Paláncz, 2016 ). These types of models have two main features: complexity and accuracy. Generally, given a certain data set, the process starts with the determination of very simple but inaccurate models. With time, more accurate but also more complex models are obtained. To prevent adverse modelling of measurement errors, data noise, or deviation, a model should be taken as a compromise between complexity and accuracy. Despite symbolic regression leading to numerous solutions, i.e. multiple equations, occasionally relatively simple equations are found. Following this approach using the software package Eureqa ( Nutonian, 2019 ), Eq. (26) is an example of such a simple equation and has been used for modelling purposes.
EUR model (Eureqa symbolic regression model) When in steady-state homogeneous fluidisation the frictional pressure drop Eq. (14) is combined with the classical drag relations or with the proposed new drag relations Eqs. (21) , ( (23) , and (25) ) in combination with Eq. (24) as well as Eq. (26) , the voidage in fluidised beds can be calculated. In the new expressions, besides the Reynolds number ( Eq. (4) ) also the Froude number is required ( Eq. (19) ).
It should be noted that all proposed models have different parameters c i .

Materials and methods
The experimental setup is presented in Section 3.1. Particle selection and fluidisation experiments are given in Section 3.2.

Experimental setup
Expansion experiments for several materials were carried out at three locations: in Waternet's Weesperkarspel drinking water pilot plant located in Amsterdam, the Netherlands; at the University of Applied Sciences Utrecht, the Netherlands; and at Queen Mary University of London, United Kingdom. In the experiments, locally produced drinking water was used. The setup ( Fig. 1 ) consisted of a 4-metre transparent PVC pipe with an inner diameter of 57 mm. Water temperature was regulated with a boiler, a cooler, and a thermostat by recirculating water through a buffer vessel connected to a water reservoir. An overflow at the top of the reactor returned water to the buffer vessel. From the buffer vessel, water was pumped through the reservoir connected to the thermostat which was set to a programmed water temperature. During terminal settling experiments, the water pump was turned off.

Particle selection and fluidisation experiments
In this study, we initially examined calcite pellets (100% CaCO 3 ) applied in drinking water softening. Polydisperse calcite pellets were sieved and separated in order to acquire more uniformly dispersed samples. To investigate hydrodynamic behaviour, also highly monodisperse and almost spherical glass beads were used: two transparent and four opaque solid-coloured SiLibeads glass beads type P. For validation purposes, liquid-solid fluidisation data were retrieved from standard references in the literature for a wide range of different particles in fluid water systems ( Wilhelm and Kwauk, 1948 ;Lewis et al., 1949 ;Richardson and Zaki, 1954 ;Loeffler, 1953 ).
Detailed information regarding particle and fluid characterisation, standard operating procedure of the fluidisation expansion experiments as well as data tabulation can be found in the Supplementary Material.

Results and discussion
Expansion experiments and fluidisation characterisation observations are given in Sections 4.1 and 4.2 . The determined drag with persisting deviation is discussed in Section 4.3 . The experimental data leading to the estimation of the Kozeny coefficient are presented in Section 4.4 . Accordingly, the drag relations based on the Reynolds number and hydraulic models based on the Reynolds and Froude numbers are discussed in Sections 4.5 and 4.6 , respectively. Finally, the voidage prediction results are given in Section 4.7 .

Expansion experiments
The acquired experimental data set consisted of a matrix with varying temperatures, grain sizes, and flow rates, as was required for a comparison of the theoretical fluidisation models. In total, 97 fluidisation experiments were carried out for calcite pellets (61) and glass beads (36) which were compared to fluidisation characteristics obtained from the literature (42). Fig. 2 shows, as an example, a typical expansion curve in which the voidage and pressure difference along the whole bed was measured for increasing superficial velocities at different temperatures. For increasing temperature, total differential pressures in the fixed bed and bed voidage in the fluid state decrease due to the decrease in viscosity i.e. less interaction force on each particle brought about by the fluid.

Fluidisation characterisation observations
During fluidisation experiments with calcite grains, we observed open spaces of water between the fluidised particles, even at relatively low fluid velocities ( Fig. 3 ). Moreover, significant heterogeneous particle-fluid patterns were detected at higher fluid velocities, even though in the literature liquid-solid fluidisation systems are generally considered to be homogeneous at lower superficial fluid velocities ( Di Felice, 1995 ). Supporting videos showing liquid-solid fluidisation experiments are shared in the Supplementary Material ( Kramer, 2019 ) to visualise the voids and bubbly flow aspects for both calcite pellets and glass beads.
To eliminate the influence of particle shape and polydispersity, the experiments were repeated with highly spherical monodisperse glass beads and we saw that the heterogeneous flow phenomena emerged yet again. The voids of water were found to provide pathways of lower resistance to the fluid, resulting in a lower drag compared to the assumed drag for homogeneous fluidisation. Figs. 4 and 5

Fig. 2.
A typical expansion experiment with SiLibeads glass beads (d p = 2.5 mm), (6 < T °C < 27), with measured differential pressure along the whole bed and voidage against the superficial fluid velocity.

Persistent drag deviation
The experimentally determined drag coefficients f T and f L for calcite pellets, glass beads, and data obtained from the literature are plotted in Figs. 6 and 7 . Classical models (Kozeny, Ergun, Carman-Kozeny and van Dijk) have been added; they are not constrained by boundary conditions to emphasize the effects of laminar and turbulent behaviour. Wall effect corrections have not been implemented because they are insignificant under our experimental conditions. A further and more detailed explanation can be found in the Supplementary Material.
When the drag f T based on the Reynolds number had been determined for calcite pellets, a certain degree of spread in the data was found ( Fig. 6 ). To eliminate the influence of particle shape and polydispersity, the experiments were repeated using highly spherical and monodisperse glass beads. Yet again a similar degree of spread in the data was found. To accentuate the spread of drag measurements in the transitional region and to eliminate the loglog effect, the standard drag f T was converted into f L using Eq. (11) .  Accordingly, the quadratic fluid velocity term in Eq. (15) is converted into a linear inversely proportional relationship between f L and velocity ( Eq. (16) ) which reduces the magnitude of drag. The presence of the spread is evident in Fig. 7 .
Data obtained from the literature was added, which shows a substantially larger spread compared to the experimentally determined drag in this work. Data from the literature, however, has a wider variance due to the use of different types of grain material ( Table 1 ).

Estimation of the Kozeny coefficient
The experimentally obtained fluidisation characteristics for glass beads in fixed bed states was used to calculate the Kozenycoefficient K . The drag relation Eq. (13) given by Erdim et al. (2015) ,   which is based on the differential pressure, can be used. In Fig. 8 , f L is plotted against Re ɛ where the original coefficients were used ( Carman, 1937 ). In the vicinity of the fixed bed (voidage close to 0.4), the deviation is substantial which can be attributed to the fixed packing irregularity. The plot shows that the original value K = 180 is slightly too high. Using non-linear curve fitting, the optimal Kozeny-coefficient K gives a value close to 150 which agrees with Ergun (1952) , Yang (2003) and many others, despite the fact  that the Kozeny-coefficient K does not have a fixed value for different conditions ( §2.1.5). Regarding spherical calcite pellets, the deviation is even higher compared to glass beads due the influence of irregularly shaped grains. To continue the modelling, a fixed value of K = 150 is assumed henceforward. Graphs and data regarding calcite pellets are given in the Supplementary Material.

Drag relations based on the Reynolds number
First, the results for classical models will be discussed. The experimentally obtained fluidisation characteristics data were used to calculate the drag ( f T and f L ) as a function of the modified particle Reynolds number Re ɛ , which is presented in Table 2 and graphs are given in the Supplementary Material. In general, when f T is used, i.e. in the laminar flow regime, drag at lower Reynolds values is For all data regarding the fluidised state, the Kozeny model violates the laminar boundary conditions and can therefore not be used to predict the drag or bed voidage. In general, the majority of liquid-solid fluidisation flow regimes occur in the transition region, in particular when the fluid is water. Exclusively laminar and turbulent flow regimes barely occur. Therefore, the Burke-Plummer model is not suitable for prediction purposes either. Conversely, the Ergun model is relatively accurate at low Reynolds numbers, but for higher Reynolds numbers, drag is increasingly overestimated. Fitting the Ergun parameters increases the prediction accuracy.
When Carman-Kozeny is examined for both calcite pellets and glass beads, drag is slightly overestimated. The most obvious explanation is that the Kozeny constant K is not constant and more likely approaches 150 rather than 180. In addition, Carman (1937) mainly used experimental data based on gas-solid systems to calibrate the Carman-Kozeny model parameters ( Fig. 9 ). A minority of Carman's dataset ( Green and Ampt, 1911 ) consisted of a few water-solid experiments at laminar conditions. Another disadvantage of this semi-empirical model lies in the assumptions that the particles are perfectly round, and that the particle distribution remains homogeneous. Since these assumptions are not fulfilled in practice, the model parameters must be adjusted, often empirically, to increase the voidage prediction accuracy. In other words, the Carman-Kozeny model parameters often used in the literature are not automatically suitable for liquid-solid fluidisation. This supports our attempts to use a fitting method to find the most reliable model parameters based on our experiments. In Table 2 , the fitted Carman-Kozeny equation has the highest R 2 . The van Dijk model ( van Dijk and Wilms, 1991 ) is valid in the transitional region but provides very diverse correlation coeffi-cients ( f T and f L ). However, fitted van Dijk parameters provide reasonable R 2 values but cannot cope with the apparent deviation. In summary, drag prediction based merely on the Reynolds number in classical models is only accurate to some extent when the parameters are fitted based on experimental data. The spread of data in terms of drag, however, is not resolved. This can also be seen in the drag-plots presented in literature ( Hoyland, 2017 ) where the data spread is less visible due to the use of logarithmic scales.
If the experimental data set is examined for glass beads using f L , the lowest K -value of 150 is confirmed. We observe that the lowest determined drag for glass beads and calcite pellets corresponds to f T = 0.86, which is lower than the minimum of 1.75 proposed in the literature ( Burke and Plummer, 1928 ;Ergun, 1952 ). The lowest determined drag regarding data obtained from the literature is even lower: f T = 0.38. These lower drag values are most likely caused by non-homogeneous fluidisation characteristics. This explains why a porous media model is less accurate and not suitable for use. Based on the numerical results in Table 2 , classical drag relations with only a Reynolds relationship are less accurate and therefore less suitable for water treatment processes.

Hydraulic models based on the Reynolds and Froude number
The persistent spread and deviations (in f T and f L ), in addition to our visible observations, reinforce our hypothesis that drag cannot be estimated accurately as a function of the Reynolds number only. The main reason for this is that crucial information about the fluidisation quality is missing. To take into account laminarturbulent as well as homogeneous-heterogeneous fluidisation characteristics, an improvement of the drag relation for liquid-solid fluidisation is proposed based on the Reynolds and Froude numbers. The effect of the extra dimension, expressed by the additional Froude number, can be seen in the 3D plot in Fig. 10 for calcite pellets. 3D plots for glass beads and data from the literature are given in the Supplementary Material.
For all three proposed types of models (cf. Section 2.2 ), K = 150 was used. In general, the correlation coefficients for glass beads, Fig. 9. Original data extracted from Carman (1937) : correlation for beds of spherical particles ( Fig. 1 ).   Fig. 10. Rotated 3D plot (f L versus Re ε , Fr p ) based on experimental data (calcite pellets) in the fluidised state. calcite pellets, and data from the literature are at least 0.95 or higher ( Table 3 ). In case the original Carman parameters are used as input parameters, in all cases, the R 2 is substantially lower. For many models, the highest obtained correlation coefficient was R 2 = 0.993, which indicates the quality of the experimental data set and model. Both RIO 1 and 2 models are based on the Carman-Kozeny equation. The RIO 1 model has three extra parameters and has the highest R 2 value, while the RIO 2 model only has one extra parameter with a slightly lower R 2 value. Although the EUR model is relatively simple and has only one extra parameter and a boundary condition, it still has a reasonable R 2 .
Regarding the RIO 2 model, the transition to turbulence is expected to occur when the two terms in Eq. (25) have approximately the same value, i.e. when the Reynolds number exceeds Re ɛ > 150 0-20 0 0 (for the obtained value of the parameter c 4 ), in accordance with Ergun (1952) .
In Fig. 11 , the drag coefficient f L is plotted against the Reynolds number through a linear plot. In this figure, the effect of the projection of the 3D plot shown in Fig. 10 on a single 2D plane, leading to the apparent spread of the drag versus modified Reynolds number, is clearly visible. The influence of the Froude number becomes apparent in Fig. 12 .   Fig. 12. Drag coefficient f L versus Fr p for glass beads and RIO 1 model (highest R 2 ).

Voidage prediction
The investigated drag models were used to predict the voidage for the fluidised state using the differential pressure Eq. (15) . The models were compared with experimental data: glass beads and calcite pellets as well as data obtained from the literature. Results are presented in Table 4 . The prediction accuracy for glass beads and calcite pellets is roughly the same. For all models, the prediction accuracy for data from the literature is lower.
Although the laminar flow regime-based Kozeny model seems to predict the voidage quite well, the model is barely valid when the given boundary conditions are respected. This is also the case with the turbulent flow regime-based Burke-Plummer model. The Ergun model is valid for all flow regimes but has a lower accuracy, i.e. a higher average relative error, compared to the original Carman-Kozeny equation. The latter model, however, is not valid for the whole regime. The van Dijk model, often applied in water treatment processes, is a little more accurate compared to Ergun. It is, however, only valid for low transition flow regimes.
During the experiments, the fluid flow was increased until voidage values were attained of ɛ ≈ 0.95. This means that the developed models are valid up to ɛ ≈ 0.95. The prediction inaccuracy is 0.8-1.5% for glass beads, 1.2-1.6% for calcite pellets, and 3.4-4.4% for data from the literature. The voidage prediction accuracy for glass beads, calcite pellets, and data obtained from the literature as given in Table 4 is based on fitting parameters based on perfectly monodisperse glass beads. For other specific grain types, besides these glass beads, it is nevertheless possible to increase the prediction accuracy in case specific fit parameters are used in the models. The most accurate model is the Carman-Kozeny-based RIO 1 model, with an overall average relative error below 1. For homogeneous fluidisation, where the Froude number has a minor effect, the models scale to the Reynolds number. For heterogeneous fluidisation, the Froude number has more influence and increasingly reduces the drag.
The choice for one of the four presented models depends on criteria such as correlation coefficient, voidage prediction accuracy, familiarity, required boundary conditions, number of fit parameters, simplicity, applicability, and/or user preferences. Despite the fact that the Carman-Kozeny-based RIO 2 model has a slightly lower voidage accuracy for glass beads and slightly higher for calcite pellets, it has two fit parameters less than the RIO 1 model and it meets all selection criteria. This means that this model is slightly more preferable.
In future work, thorough CFD modelling is recommended to determine whether the Froude number is a suitable dimensionless number to exactly describe or evaluate the concept of 'heterogeneous flow phenomena'.

Conclusions
Liquid-solid fluidisation processes are frequently used in industry, such as drinking water treatment processes. For pellet softening, the operational field lies in the vicinity of incipient fluidisation to provide a large crystallisation surface area and consequently to obtain optimal process conditions. This operational field falls within the initial transitional flow regime rather than turbulent flow regimes. To obtain optimal process conditions, the overall fluidised bed voidage is a crucial process parameter which can be estimated by means of drag relations. Traditionally, in drag relations, the emphasis of the dimensionless standard drag coefficient, defined as a function of the Reynolds number, is focused on the turbulent flow regime ( ~v s 2 ) and less so on the laminar ( ~v s ) and transition flow regimes.
We propose four adjustments to improve the drag analysis to mathematically describe the fluidisation stability and to increase the overall voidage prediction accuracy. The first and second adjustments are to multiply the standard drag coefficient with the Reynolds number and to use a double linear representation instead of a traditional double logarithmic representation; this improves the distinctive capability of the drag analysis considerably. The third adjustment concerns coping with heterogeneity phenomena in liquid-solid fluidised beds. The traditional drag relation based on the particle Reynolds number is extended with the particle Froude number. By adding a third dimension to the traditional 2D linear plots, where drag is plotted against Reynolds, the apparent and persistent spread in effective full-scale drag can be explained and visualised by means of a 3D plot. The fourth adjustment increases the accuracy of the dimensionless drag coefficient and therefore also the voidage prediction by using experimental data based on total mass of particles, bed height, and particle density measurements rather than using sensitive differential pressure measurements.
Four new prediction models have been synthesized (SON, RIO 1, RIO 2, and EUR) which enable us to predict the voidage in the fluidised state more accurately. The prediction average relative error decreased from approximately 5% using the best literature equation (exclusively based on Reynolds number) to 1-2% with the new equation (based on Reynolds and Froude numbers). The RIO 2 model based on Carman-Kozeny has a voidage prediction inaccuracy of only 1% and can be used for calcite pellets as well as for spherical grains used in full-scale drinking water treatment processes such as pellet softening.
Supplemental Material: Supplemental data for this article can be accessed at doi…

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.