Investigation of drag models in CFD modeling and comparison to experiments of liquid–solid fluidized systems

International Journal of Mineral Processing 104-105 (2012) 58–70

Contents lists available at SciVerse ScienceDirect

International Journal of Mineral Processing journal homepage: www.elsevier.com/locate/ijminpro

Investigation of drag models in CFD modeling and comparison to experiments of liquid–solid fluidized systems Olli Visuri a,⁎, Gijsbert A. Wierink b, Ville Alopaeus a a b

Aalto University, School of Chemical Technology, Department of Biological and Chemical Engineering, PO Box 16100 FI-00076 Aalto, Finland Aalto University, School of Chemical and Technology, Department of Materials Science and Engineering, PO Box 16100 FI-00076 Aalto, Finland

a r t i c l e

i n f o

Article history: Received 23 May 2011 Received in revised form 22 December 2011 Accepted 27 December 2011 Available online 5 January 2012 Keywords: Euler–Euler Drag Fluidization CFD Digital imaging

a b s t r a c t A liquid–solid fluidization system was investigated with Computational Fluid Dynamics (CFD) by using a transient Eulerian-Eulerian model. The study focused on various drag models between the phases and how the results vary when simulating the system 2-dimensional or 3-dimensional. Also the grid dependencies to the results were investigated. The simulation results were validated experimentally using a digital imaging method. The suitability of this experimental method is also investigated in this paper. The CFD results show that the outcome from different drag models vary considerably and therefore the used model has to be chosen with care. © 2012 Elsevier B.V. All rights reserved.

1. Introduction Multiphase reactors are very often used in chemical industry. These involve gas–liquid, liquid–liquid, liquid–solid and gas–liquid– solid phases. In order to optimize the reactor design in industrial scale it is important to understand the hydrodynamics of the reactor. Traditionally this has been done with experimental scale-up. However, the correlations gained this way depend on the vessel geometry and therefore the final geometry has to be specified fairly early during the process development. If several options on the vessel dimensions are investigated the increase of costs is usually high. In recent years due to the increase of computational capacity, Computational Fluid Dynamics (CFD) has been applied extensively to model these reactors and their flow behavior. CFD in multiphase systems can be divided into Eulerian–Eulerian (E–E) or Eulerian– Lagrangian (E–L) modeling. In the E–E method both phases are described using the Eulerian conservation equations and the phases are assumed to be interpenetrating continua (Drew, 1983), (Bothe et al., 2007). These are represented as averaged conservation equations. In E–L modeling the carrier phase is still described in an Eulerian framework while the dispersed phase is described using a Lagrangian frame of reference. For CFD simulation of more dense systems where concentration profiles are a key result, the E–E method is preferred (Enwald et al., 1996).

⁎ Corresponding author. Tel.: + 358 470 25813. E-mail address: Olli.Visuri@tkk.fi (O. Visuri). 0301-7516/$ – see front matter © 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.minpro.2011.12.006

Although CFD is a promising tool for designing the equipment, there is a need to validate the CFD results, since the phenomena inside the cell volume has to be modeled. The goal in CFD is to be independent on the system, but the closure models that are used are commonly direct derivations from certain experimental system. Therefore their consistency should be checked experimentally before using them in modeling industrial systems. Normally the models are validated with much simpler systems, using so called mock-up experiments, where the chemicals are safer and therefore the need to invest complex systems during laboratory and pilot scale is reduced. For experimental validation for multiphase systems numerous different methods can be found from literature. Basically these can be divided into two categories, invasive and non-invasive. In invasive techniques, commonly a probe (e.g. conductivity, thermal, optical) is placed inside the system. These methods provide local data from the system, especially when considering spatial and time resolution. The drawback with these systems is that they tend to disturb the experiment. Also gaining data from multiple locations causes problems. In non-invasive techniques the analysis is made outside the system, commonly using visualization or imaging technique. The drawbacks with these techniques are usually that they are limited to low holdups of the dispersed phase and the need of transparent equipment. A comprehensive review of various experimental techniques is presented (Boyer et al., 2002). In this work the main focus was to investigate how a liquid–solid fluidized bed can be modeled with CFD. The experimental system consists of rectangular pipe where the solid bed height and concentration is measured using photographic technique. The fluid flow was driven through a dense grid making the turbulence isotropic.

O. Visuri et al. / International Journal of Mineral Processing 104-105 (2012) 58–70

Total amount of particles placed into the system was measured and the particle distribution was analyzed with photographic method. This method is non-intrusive, so there is no disturbance to the flow field. The numerical simulation was carried out using the open source CFD package OpenFOAM, version 1.6, released by OpenCFD Ltd. (OpenCFD Ltd., 2009). The system was modeled by OpenFOAM's Eulerian-Eulerian solver twoPhaseEulerFoam, based on the work of Rusche (2002). 2. Experimental 2.1. Experimental setup and procedure 2.1.1. Particles The used particles were porous cylindrical extrudates (see Fig. 1a), which are common in fluidized beds in chemical industry. The porosity of the used particles was 0.4744. The density and size distribution was measured by photographing the particles. It was assumed that the cross sectional area is the same for every particle and the length of the particles were determined from the image. The wetted density

of the particles was approximately 2173 kg/m 3, mean volume of the surf aceareaequivalent particles were 1.8 mm 3 and the average sphericity ( surf acearea ) real were approximately 0.81. The volumetric size distribution is presented in Fig. 1b. 2.1.2. Experimental apparatus The experimental apparatus consisted of a transparent vertical channel of square cross-section that was used for the fluidization. The particles mentioned above were suspended using water. A schematic diagram of the apparatus is shown in Fig. 2. The dimensions of the pipe were 0.1 m (width) 1.8 m (height) and 0.1 m (depth). 0.6 m from the bottom a dense metal grid (internodal distance approx. 0.7 mm) was placed and the particles were placed above this grid. A 2 mega pixel digital video camera (JAI, CV-M2) with 0.95/ 25 mm lens (Schneider optics, Xenon) was used in the experiments. At the acquisition the f-number was 2.0. Shutter speed was set to 1/ 250 s. The camera's frame rate was 30 frames/s. The system was illuminated with a 2 kW HQI lamp using it as backlight. A set of diffusers were used to equalize the light intensity over the whole tube.

a

b

800 700

No. of Particles

600 500 400 300 200 100 0 0

1

59

2

3

4

V [mm3] Fig. 1. An image of the particles and the volumetric distribution of the particles.

60

O. Visuri et al. / International Journal of Mineral Processing 104-105 (2012) 58–70

Fluid outlet

z x

Camera

y

Particles

Grid

A 2 kW HQI lamp A set of diffuser films Fluid inlet Fig. 2. A schematic view of the experimental setup.

2.1.3. Experimental procedure The experiments were done with three different particle amounts and each one of them was fluidized with four different velocities, totaling 12 different experimental conditions (see Table 1). From each condition a set 150 images were acquired and the experiments were repeated 3 times and the average values of these were used. For the bed height determination the experimental error was estimated to be approximately ±2 cm. This was analyzed manually from the graphs produced from the experiments. 2.2. Digital imaging method A digital image can be also considered as a numerical matrix. Element of the matrix is a pixel that describes the light intensity of a point of the image (in color images the element is a sub matrix). A typical grayscale image uses 8-bits to describe the intensity of the pixel, thus the value can vary from 0 to 255 (equals 2 8), where 0 equals to black and 255 to white. When using backlight as we did in our experimental apparatus, the particles will show as dark spots in the image. The images were processed with a java-based ImageJ-freeware (ImageJ-software, Research Services Branch, National Institute of Health, USA). Additional user-implemented analysis routines were developed for the image analysis. 2.2.1. Background removal Before calibration and experimental analysis the imaging background of the system was established by taking 5 images running the system without particles. These images were analyzed and the results were averaged to obtain the background values. The

Table 1 Experimental conditions. Mass of the particles (g)

Flow velocity (m/s)

396.2

0.031 0.058 0.085 0.100 0.032 0.061 0.080 0.091 0.033 0.058 0.078 0.093

635.0

791.7

background was removed from the actual experiments during calibration and data analysis. 2.2.2. Calibration The calibration was done by running the apparatus with a predetermined amount of particles and fluidizing it as uniformly as possible. Images of the calibration were captured with the camera and the average intensity was calculated. The intensity value was averaged over the whole target. Naturally when the particle volume was increased at the tube some inaccuracy occurred, since the bed did not distribute evenly into the system. For the actual calibration between the light intensity (I) and concentration (c) it is possible to apply Lambert–Beer law (Atkins (1998)) I

∫ I0

b

dI ¼ ∫ −ac dz I 0

ð1Þ

where a is a coefficient, I the intensity, b height of the measured system and z represents the direction of height. Originally the Lambert–Beer law is applied in spectroscopy and the basis is in the reduction in the intensity that occurs when light passes through a layer of molecules, which are selectively absorbing certain wavelengths. The analogy in this experiment is similar, since the orientation of the particles is random in the system. Also the averaged volume is large enough that the particles do not block the whole measured area within the feasible measuring limit. By replacing I with Imeasured − Iref and I0 with Imax − Iref the measuring limits that comes from the imaging analysis are set to the calibration. Iref represents the minimum intensity (e.g. infinite concentration) that can be detected with this method and Imax is the maximum intensity that can be obtained from the analysis (e.g. the zero concentration). In 8-bit images it is 255 (=2 8, the range varies from 0 to 255). Intergration of Eq. (1) with these end points yields Eq. (2) c¼−

! Imeasured −I ref 1 ln : ab I max −Iref

ð2Þ

From the calibration results the values for Iref and a were fitted for each calibration. The reference value has a high impact on the standard deviation of the fitted curve and therefore it has to be chosen with care. The calibration curve is shown in Fig. 3. In Fig. 4 the deviation between the experimental results and the fitted curve is shown. From this graph it can be concluded that the experimental errors are more random than systematic. Besides analyzing the

O. Visuri et al. / International Journal of Mineral Processing 104-105 (2012) 58–70

90

100 Measured Fitted

80

1 Intensity Unit 2 Intensity Units 3 Intensity Units 4 Intensity Units 5 Intensity Units 6 Intensity Units 7 Intensity Units 8 Intensity Units

80

70 60

+20% deviation

60

50

C (g/l)

C (g/l)

61

40

40 30 20

20

no deviation

10 0 0

50

100

150

200

250

0

0

10

20

l (-) Fig. 3. Calibration between concentration and the intensity value from the images.

standard deviation the effect of random error should be checked since the calibration is highly non-linear, and may cause large errors with small changes in the intensity value. This can be done with a sensitivity analysis, where it is checked how much the variation affects on the output of the calibration correlation. In Fig. 5 is shown how much the concentration changes compared to the original value if there is error in the intensity value. The investigation was done with intensity values from one to eight. From the graph it can be seen, that the deviation starts to be exceed 20% when the concentration is approximately 35 g/l. Therefore the values that were above 35 g/l were excluded from the analysis. This limitation does not have any effects, when analyzing the bed height, and also it is possible to analyze the concentration profiles on the top of the bed.

20

40

50

60

C (g/l) Fig. 5. Sensitivity analysis for the experiments.

3. Modeling To model the fluidized bed, the open source CFD OpenFOAM-1.6 (OpenCFD Ltd., 2009) was used. OpenFOAM's meshing utility blockMesh was used to generate structured orthogonal meshes. The physical dimensions of the computational domain are 0.10 m wide, 0.10 m deep and 1.20 m high. The field for dispersed phase α was set to the same level as the experimental results gave. Fig. 7 shows an example mesh in 2D. Different mesh densities in both 2D and 3D were tested in this work. In Table 2 the mesh densities are summarized. It was assumed that the grid levels out the jet effect caused by the inlet flow of water and the inlet velocity in the simulations were set as constant. 3.1. TwoPhaseEulerFoam solver

2.2.3. The analysis The recorded video sequence was divided into individual images and further into smaller blocks which height were 15 pixels (equaling approx. 1.1 cm, see Fig. 6). The average intensity of these blocks was calculated, converted to concentration with the calibration curve and furthermore converted to void fractions as seen in Figs. 13 and 14. The parts of the image where the rack of the experimental setup was were excluded from the analysis.

The twoPhaseEulerFoam solver is an Eulerian–Eulerian solver readily available in the OpenFOAM CFD package (OpenCFD Ltd., 2009) and based on the work of Rusche (2002). Cell-centered and phase averaged equations are solved using the Finite Volume Method

10 8 6

Deviation (%)

4 2 0 -2 -4 -6 -8 -10

C (g/l) Fig. 4. Deviations between experiments and the correlation.

Fig. 6. A visualization from the block division.

62

O. Visuri et al. / International Journal of Mineral Processing 104-105 (2012) 58–70

(FVM). For the dispersed phase α, the continuity and momentum equations are (Silva et al., 2008) ∂ðαρα Þ þ ∇⋅ðαρα U α Þ ¼ 0 ∂t

ð3Þ

  ∂ðαρα Uα Þ eff þ ∇⋅ðαρα Uα Uα Þ ¼ ∇⋅ αTα þ Mα þ αρα g ∂t

3.2. Drag models ð4Þ

Teff α

where the effective stress tensor for phase α is composed of the normal and deviatoric components as   2 2 eff eff Tα ¼ −pα I þ ν α 2Dα − ð∇⋅Uα ÞI − ρα kα I: 3 3

ð5Þ

In Eq. (5) all parameters with index α refer to phase α, I is the identity matrix, D is the strain rate tensor, ν is the kinematic viscosity, k is the turbulent kinetic energy, p the pressure and ρ the density. The momentum source terms in Eq. (4) are the divergence of the interfacial stress tensor, a momentum coupling term M and a body force term. The momentum source term M in Eq. (4) was derived in its averaged form by Gosman et al. (1992) and contains a drag coupling term Mα,drag of the form (Rusche, 2002)  ð1−α  ÞKUr Mα;drag ¼ α

where Ur is the relative velocity. Depending on the formulation of the (stationary) drag coefficient CD, Ur is either the apparent relative velocity αg(Ug–Up) (Boemer et al., 1995) or the relative interstitial velocity Ug–Up (Enwald et al., 1996).

ð6Þ

In literature there are hundreds of different drag models that have been introduced since Stokes has introduced the modeling. Drag models are commonly based either on the pressure drop over a fluidized bed or on the drag coefficient on a single particle (Enwald et al., 1996). The drag function K is then calculated based on the appropriate model for pressure drop or drag coefficient on a particle. In this work the idea was to collect models that have different approach to the drag force. The models used here are models that have been commonly used for decades, such as the models by Ergun (1952), Schiller and Naumann (1933), Wen and Yu (1966) and Syamlal and O'Brien (1987). Also two combination models were checked, developed by Gidaspow (1994), where the model changes depending on the conditions inside the cell volume and two models which take into account the non-sphericity of the particles (Haider and Levenspiel, 1989; Swamee and Ohja, 1991). Although initially the tested models have been developed in a specific system and conditions the usages of these models have been extended outside their initial conditions (e.g. Reddy and Joshi (2009) Asif and Ibrahim (2002), Panneerselvam et al. (2007)). 3.2.1. Schiller–Naumann Stokes developed an analytical expression for the stationary drag coefficient CDs on a sphere for Reynolds number (Re) ≪ 1 as C Ds ¼

24 : Re

ð7Þ

By including inertia, Oseen (1913) formulated an expression that extends Eq. (7) to flow regimes characterized by a higher Reynolds number as C Ds ¼

  24 3 1þ Re : Re 16

ð8Þ

Eq. (8) is in fact a truncated form of the Goldstein expansion (Goldstein, 1929). Schiller and Naumann (1933) found that for higher Reynolds numbers the Oseen equation Eq. (8) deviates from experimental observations for falling spheres. For the regime of Re b 800, Schiller and Naumann (1933) wrote the drag coefficient (CD) as CD ¼

 24  0:687 1 þ 0:15Re : Re

ð9Þ

The Schiller–Naumann equation, Eq. (9), predicts a lower drag coefficient than the Oseen equation, Eq. (8). For Re ≈ 1000 Eq. (9)

Table 2 Mesh densities with physical dimensions x = 0.10 m, y = 0.10 m, and z = 1.20 m, with total number of computational cells Ncells. Dimension

Fig. 7. Example of a 2D mesh of the computational domain as generated using OpenFOAM's blockMesh utility. The void fraction of dispersed phase is indicated by “alpha”.

2D 3D 2D 3D 2D 3D 2D

Cartesian axis

Ncells

x

y

z

10 10 20 20 40 40 100

1 10 1 20 1 40 1

120 120 240 240 480 480 1200

1200 12,000 4800 96,000 19,200 768,000 120,000

O. Visuri et al. / International Journal of Mineral Processing 104-105 (2012) 58–70

coincides with the drag coefficient of 0.44 for a smooth sphere. Therefore, the Schiller–Naumann drag coefficient becomes ( C Ds ¼

 24  0:687 1 þ 0:15Re ; Re≤1000 Re 0:44; Re > 1000

ð10Þ

and the drag function is K¼

3C Ds ρl U r 4dp

ð11Þ

F g −F s F ks

ð12Þ

where Fg is gravitational force, Fs is buoyancy force of fluid acting on the particle and Fks is the drag force of the fluid acting to the particle. The investigation was based on experiments done with various liquid–solid systems. The velocity ranged from laminar to turbulent conditions and they compared their results to four other data sets. Wen and Yu (1966) also studied the minimum fluidization velocity with various gas–solid and liquid–solid systems and compared their results to 16 different data sets. Using a correlation to fit data of Richardson and Zaki (1954), Wen and Yu (1966) wrote a drag function as K¼

   −2:65 3 C α 1−α p ρl U r 1−α p 4dp Ds p

ð13Þ

where CDs is the Schiller–Naumann drag coefficient, Eq. (10). 3.2.3. Ergun Ergun (1952) devised a combined equation valid for both low and high Reynolds number flows in a packed bed. Ergun (1952) considered a viscous regime for flows characterized by a low Reynolds number, where the Kozeny–Carman (Kozeny, 1927; Carman, 1937) equation shows good agreement with experimental data. For higher Reynolds numbers, however, the kinetic Burke–Plummer (Burke and Plummer, 1928) equation shows better agreement. The form of the Ergun equation is based on the loss of energy over the packed bed, due to friction. Here, Ergun used the Blake friction factor of the form (Ergun, 1952) f ¼ 150

αp þ 1:75: Re

ð14Þ

Ergun (1952) assumed that all energy dissipates according to the friction factor Eq. (14) and wrote the ratio of total loss of energy to the loss of kinetic energy as  Δpgdp

1−α p

Lρl U 2r

αp

and replacing the pressure drop per unit length by the pressure gradient, the Ergun equation becomes (Enwald et al., 1996) K ¼ 150

3.2.2. Wen–Yu The model created by Wen and Yu (1966) was originally validated based on two types of experiments. Wen and Yu (1966) studied the voidage function f ðεÞ ¼

In Eq. (16) the left term on the right hand side represents viscous part of pressure loss, i.e. the Kozeny–Carman equation, and the right term on the right hand side represents the kinetic component, i.e. the Burke–Plummer equation for turbulent flows (Ergun, 1952; Gibilaro et al., 1985). By writing the drag function as   1−α p ∇p K¼ ; ð17Þ Ur

2

where dp represents the particle diameter.

ð15Þ

where L is the characteristic length of the fluidized bed. Equating the friction in Eq. (14) with the relative kinetic energy in Eq. (15) Ergun obtained α 2p αp Δp μU ρ U2 ¼ 150 l r   þ 1:75 l r   : L dp g 1−α 3 dp g 1−α 3 p

p

ð16Þ

2

αp αp μl ρU   þ 1:75 l r   : dp 1−α 2 dp 1−α 2 p

ð18Þ

p

Ergun's work is based on gas–solid systems in fixed beds. The main focuses for determination were the rate of the fluid flow, viscosity and density of the fluid, closeness and orientation of the solids and the size, shape and surface of the solid particles. 3.2.4. Gidaspow–Ergun–Wen–Yu The Wen–Yu model was developed for particulate systems where viscous forces are dominant, i.e. relatively dilute flows (Lundberg and Halvorsen, 2008). Therefore, Gidaspow combined the Wen–Yu and Ergun drag models for relatively dilute and dense flows, respectively. The Gidaspow–Ergun–Wen–Yu model uses the Ergun equation, Eq. (18), for αp ≥ 0.2 and the Wen–Yu equation, Eq. (13) for αp b 0.2. van Wachem (2000), however, points out that the step-change at αp = 0.2 can lead to numerical instability. 3.2.5. Gidaspow–Schiller–Naumann The Schiller–Naumann drag model, Eqs. (10)–(11), uses a drag coefficient for a single spherical particle. To account for particle drag in a multiparticle system, Gidaspow and Ettehadieh (1983) used the voidage function of the work of Richardson and Zaki (1954). The Gidaspow–Schiller–Naumann drag model is (Gidaspow and Ettehadieh, 1983; Enwald et al., 1996): K¼

   −2:65 3 C α 1−α p ρl U r 1−α p ; 4dp Ds p

ð19Þ

where CDs is as defined in Eq. (10), with the modification that Re in Eq. (10) is replaced by (1 − αp)Re (Rowe, 1961). Eq. (19) is however only valid for dilute systems with αp ≤ 0.2 (Enwald et al., 1996). 3.2.6. Gibilaro For higher Reynolds numbers Gibilaro et al. (1985) found that the Ergun equation, Eq. (17), significantly deviates from experimental data for friction factor in a packed bed. MacDonald et al. (1979) and Gibilaro et al. (1985) write that the Blake–Kozeny constant, 150 in Eq. (18), should be 180 to fit experimental data better and that the Ergun equation is not valid for the limit of a fully expanded bed. To account for very low values of αp and higher Reynolds numbers, Gibilaro et al. (1985) rewrote Eqs. (17) and (18), giving (Enwald et al., 1996)  K¼

3

63

 −2:8 17:3 ρ jU j  þ 0:336 l r α p 1−α p : Re dp

ð20Þ

The model was validated with pressure drop experiments using either gas or liquid as the fluid. The experiments were done in fixed bed with either a cubical array or random packing. The voidage varied from 0.41 to 0.94. 3.2.7. Syamlal and O'Brien The Syamlal and O'Brien (1987) drag model is based on the notion that for single- and multiparticle systems the Archimedes number remains the same for terminal settling velocity (Ramesh and

64

O. Visuri et al. / International Journal of Mineral Processing 104-105 (2012) 58–70

Raajenthiren, 2010). Syamlal and O'Brien (1987) use a simplified form of an expression of Garside and Al-Dibouni (1977) for the relation between settling velocity and void fraction V r −A ¼ 0:06Res ; B−V r

ð21Þ

where Res is the Reynolds number for a single particle. Vr is the ratio of terminal settling velocity in a multi-particle system to that of a single particle system. Syamlal and O'Brien (1987) then substitute Res = Re/ Vr in Eq. (21) and obtain  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi V r ¼ 0:5 A−0:06Re þ ð0:06ReÞ2 þ 0:12Reð2B−AÞ þ A2  4:14 A ¼ 1−α p 8  1:28 > < 0:8 1−α p ; α p ≥0:15 2:65 : B¼  > : 1−α p ; α p b0:15

Fig. 8. The principle axis of a particle for the Corey shape factor.

ð22Þ

where β is the Corey shape factor (see Fig. 8) c β ¼ pffiffiffiffiffiffi ab

Using the drag coefficient for a single particle of Dallavalle (1948) !2 4:8 ð23Þ C Ds ¼ 0:63 þ pffiffiffiffiffiffiffiffi Res Again substituting Res = Re/Vr, the drag function in the Syamlal– O'Brien drag model (Benyahia et al., 2005) is:   sffiffiffiffiffiffiffiffiffiffiffiffi  ffi!2 3α p 1−α p ρg Vr K¼ 0:63 þ 4:8 ð24Þ jUr j: Re 4V 2r dp 3.2.8. Haider–Levenspiel Most of the models mentioned above have been developed for spherical particles. Haider and Levenspiel (1989) developed a model which takes into account the sphericity (ϕ) of the particle. Sphericity is defined as the ratio between volume equivalent surface against the actual surface of the particle. rffiffiffiffiffiffi 3 9π 2 =3 φ¼ ð25Þ V A 16

ð29Þ

where a b b b c. This model has been validated using two groups of unstructured particles, crushed rock and natural sediments. Since the experimental work has been done with β values of 0.3, 0.5, 0.7 and 1.0, the model is also stated applicable when 0.3 b β b 1. This model was also implemented to the general drag function (11). 3.3. Post processing A concentration profile similar to the experimental data was extracted from the CFD results, to make comparison possible. To generate this profile from both 2-D and 3-D simulations, the volume weighted cell values for solid phase fraction were summed according to: α p ðzÞ ¼

∑y ∑x α p ðx; y; zÞdxdy ∑y ∑x V cell dxdy

;

ð30Þ

where αp(x,y,z) is the solid void fraction field, αp(z) is the solid void fraction field as function of only bed height, and Vcell is the computational cell volume.

The actual coefficient is presented as C D;0 ¼

 24  C B 1 þ A⋅Rep þ ; Rep 1 þ ReD

4. Results and discussion ð26Þ 4.1. Experimental results

p

where,   A ¼ exp 2:3288−6:4581ϕ þ 2:4486ϕ2 B ¼ 0:0964 þ 0:5565ϕ   C ¼ exp 4:905−13:8944ϕ þ 18:42222ϕ2 −10:2599ϕ3   2 3 D ¼ exp 1:4681 þ 12:2584ϕ−20:7322ϕ þ 15:8855ϕ :

ð27Þ

mi;corr ¼

This drag coefficient was implemented into the general drag function (Eq. (11)). The parameters A, B, C and D are polynomial fits from experimental data. They have been minimized by using least-squares fit. The experimental data contained totally 419 data points with 4 different sphericities. 3.2.9. Swamee–Ohja Swamee and Ohja (1991) used another approach to describe the sphericity. 48:5

0:8 ⋅Rep 0:64 1 þ 4:5⋅β0:35 " # !0:32 Rep 1 þ ⋅ 18 Rep þ 100 þ 100β β þ 1:05⋅β0:8

4.1.1. Mass balance matching To receive comparable results and due to the non-linearity of the calibration the mass (m) balance was matched with the known total mass of the particles. This was done by

C D;0 ¼

ð31Þ

Table 3 Correction factors from mass balance matching. Volume of the particles (ml)

Flow velocity (m/s)

Correction parameter

500

0.031 0.058 0.085 0.100 0.032 0.061 0.080 0.091 0.033 0.058 0.078 0.093

0.27 0.13 0.14 0.09 0.26 0.17 0.20 0.16 0.22 0.25 0.31 0.23

800

1000

ð28Þ

mreal m: ∑mi i

O. Visuri et al. / International Journal of Mineral Processing 104-105 (2012) 58–70 mreal The corrections factors (i.e.∑m ) are presented in Table 3. i As seen from the data, the factors are diverging clearly from one. This shows the limitation of the experimental method as mentioned in Section 2.2.2. In this case the bed concentration mostly exceeds the feasible concentration limit and therefore it is more comfortable to have the main focus on the validation of the bed height.

65

4.2. Comparison of CFD models and validation 4.2.1. Bed height Because of the momentary fluctuation on the top of the bed, the bed height needs to be defined on a statistical basis. The criterion was set that the bed consisted 99% of the total amount of the particles. This was determined by checking the experimental data with various

Fig. 9. Visual results from the 2-dimensional simulations, grid size 100 × 1 × 1200 cells, experimental results on the left.

66

O. Visuri et al. / International Journal of Mineral Processing 104-105 (2012) 58–70

criteria. When the criterion was above 99%, the error caused by single particles started to significantly affect the results. This indicates that the experimental error started to affect significantly on the results. Visual comparison shows that in the 2-D simulations the particles tend to concentrate on the sides of the bed. In 3-D simulations a concentrated spot is in the middle of the bed (see Fig. 9 and

Fig. 10), when simulating in 3-D and non-spherical particles with 2-D stability issues was more important. The stability could be increased by using under-relaxation during the solving process. However with time-dependant simulations the under-relaxation can have significant effect on the results, and they would not be comparable. By setting the grid coarser the convergence increased

Fig. 10. Visual results from the 3-dimensional simulations, grid size 10 × 10 × 120 cells, experimental results on the left (see Supplementary data).

O. Visuri et al. / International Journal of Mineral Processing 104-105 (2012) 58–70

a

b

0.6

0.5

0.6

0.5

0.4

0.4

h (m)

h (m)

67

0.3

0.3

0.2 0.2 0.1 0.1 0 0

5

10

15

20

0

t (s)

0

5

10

15

20

t (s)

Ergun Gibilaro Gidaspow-Ergun-Wen-Yu Gidaspow-Schiller-Naumann Schiller-Naumann Syamlal-O’Brien Wen-Yu Experiment

Haider-Levenspiel 0.81 Haider-Levenspiel 1.0 Swamee-Olja 0.3 Swamee-Olja 0.7 Swamee-Olja 1.0 Experiment

Fig. 11. The bed height as a function of simulation time. 2-D. grid size 100 × 1 × 1200 with spherical particles (a) and non-spherical particles (b). Setup 800 ml particles, liquid velocity 0.061 m/s.

significantly and the results are indicative although they might have some grid dependency. The Haider–Levenspiel model has been originally developed using spheres, cubes, tetrahedrons and octahedrons. All these shapes are such that the ratio of a sphere fitted outside the particle against the sphere fitted inside the particle is low. With a cylindrical particle this ratio is high, which gives some doubts to use this model in these simulations. Also the relationship with the sphericity and the drag coefficient is non-linear, which implies that instead of

a

considering average values, a distribution to describe the sphericity should be implemented into the model. This applies also to Swamee and Ohja's drag model. Since the modeling with twoPhaseEulerFoam is carried out in a transient mode, it is possible to follow the development of the bed height as a function of time, as seen in Fig. 11. and Fig. 12. This shows how the models start to diverge from the experimental results and shows when the bed has reached a stable point and there is no notable change in the bed height. By looking at the figures, it can be

b

0.6

0.5

0.5

0.4

0.4

0.3

h (m)

h (m)

0.6

0.3

0.2 0.2 0.1 0.1 0 0

5

10

15

t (s) Ergun Gibilaro Gidaspow-Ergun-Wen-Yu Gidaspow-Schiller-Naumann Schiller-Naumann Syamlal-O’Brien Wen-Yu Experiment

20 0 0

5

10

15

20

t (s) Haider-Levenspiel 0.81 Haider-Levenspiel 1.0 Chhabra 0.3 Chhabra 0.7 Chhabra 1.0 Experiment

Fig. 12. The bed height as a function of simulation time, 3-D, grid size 10 × 10 × 120 with spherical particles (a) with non-spherical particles (b). Setup 800 ml particles, liquid velocity 0.061 m/s.

68

O. Visuri et al. / International Journal of Mineral Processing 104-105 (2012) 58–70

Table 4 The bed heights from simulations. Model

Bed height [m] in 2D simulations

Ergun Gibilaro Gidaspow, Ergun and Wen and Yu Gidaspow, Schiller and Naumann Schiller and Naumann Syamal and O'Brien Wen and Yu Haider and Levenspiel 1.0 Swamee and Olja 0.3 Swamee and Olja 0.7 Swamee and Olja 1.0 Experimental

0.445 0.245 0.395 0.325 0.155 0.405 0.295 0.165 0.265 0.165 0.155 0.270

Bed height [m] in 3D simulations 0.275 0.412 0.375 0.202 0.335 0.205 0.385 0.230 0.195 0.270

seen that the Ergun model and model by Syamal and O'Brien did not reach stable bed height after 20 s of simulated time in the 2-D simulations. In 3-D cases these models turned the bed from fluidized bed to a moving bed, which can be seen since the bed height starts to fluctuate. This indicates that these two models predict too high velocities for particles related to the liquid in the system. For the non-spherical particles only Swamee and Ohja's model converged. With the Haider and Levenspiel model the particles start to pack on the top of the system and since boundary conditions are such that the particles cannot exit the system the simulation interrupts. This happens approximately after 7.5 s of simulated time. The bed heights after 20 s of simulation time in Table 4 are shown. It shows clearly that the tested models vary and shows that the models depend significantly on the system that they were validated with. Also a dependency of the third dimension can be noted, since all bed heights are larger in 3D-simulations, especially when considering non-spherical particles. In our validation case the model by Gibilaro shows the best match with experimental results. Those cases where the bed was strongly time-dependent (e.g. the fluidized bed turned into a moving bed) were excluded from the table.

4.2.2. Concentration profiles In Fig. 13 the predicted concentration profiles as a function of height are shown with 2-D and 3-D simulations. From the figures it can be seen that relatively low grid size gives sufficient results when estimating the bed height. In the actual concentration profiles there are differences when the grid is denser. The denser simulations show that the concentration is more uniform inside the bed. This can be noted in both 2-D and 3-D cases. Unfortunately it was not possible to gain reliable experimental data from the profile since inside the bed the concentration exceeded the limit as mentioned in Section 2.2.2. In Fig. 14 the concentration profiles from the modeling after 20 s of simulation time and the experimental points that were considered feasible as mentioned in Section 2.2.2. are shown. The modeling setup was in these cases 2-D 100 × 1 × 1200. From the figures it can be seen that the Ergun model and the model by Syamal and O'Brien over predicts the bed height and therefore the concentration profile is more likely to be too uniform. This was also noted when simulating the 3-D cases, where these two models transformed the bed into a moving bed instead of fluidized bed. The Gidaspow–Ergun–Wen–Yu model seems to have similar predictions, although no transform to moving bed occurred. This is probably due to the fact that the model uses both Ergun and Wen and Yu drag functions. The Schiller–Naumann model under predicts the bed height and in all cases the bed is predicted to collapse to the bottom of the vessel. In Gidaspow–Schiller–Naumann-model and Wen–Yu model the differences are quite small, which is probably caused by the fact, that the bed is practically modeled with Wen and Yu drag function. The model by Gibilaro predicts the results most reliably, although there is some numerical instability at the top region of the bed. 5. Conclusions In this work a digital imaging method to validate CFD simulations is presented. The method is simple and fast to validate the simulations. Also agreement between the experiments and simulations was found when monitoring the bed height. Actual concentration validation deep inside the bed was not reliable in this case, due to the non-linearity of the correlation between concentration and the intensity value of the image. Although most parts of the bed

Fig. 13. Grid dependency on the concentration profile.

O. Visuri et al. / International Journal of Mineral Processing 104-105 (2012) 58–70

0.4

0.6

0.2

0.4

0.6

0

0.4

0.6

0.4

0.6

1.2 1 0.8 0.6 0.4 0.2 0

1000 ml, 0.058m/s

0

0.2

0.4

0.4

0.6

0

0

0.6

0.2

0.4

1.2 1 0.8 0.6 0.4 0.2 0

Average void fraction (-)

0.2

0.4

0.4

0.6

800 ml, 0.091m/s

0.6

1.2 1 0.8 0.6 0.4 0.2 0

0

0.6

Average void fraction (-)

0.2

0.4

0.6

Average void fraction (-)

1000 ml, 0.078m/s

0

0.2

Average void fraction (-)

Average void fraction (-)

Height (m)

Height (m)

Height (m) 0.2

Average void fraction (-)

0.2

0.2

800 ml, 0.080 m/s 1.2 1 0.8 0.6 0.4 0.2 0

Average void fraction (-)

1000 ml, 0.033 m/s

0

0

1.2 1 0.8 0.6 0.4 0.2 0

Average void fraction (-)

Height (m)

Height (m)

1.2 1 0.8 0.6 0.4 0.2 0

Average void fraction (-) 1.2 1 0.8 0.6 0.4 0.2 0

0.6

800 ml, 0.061 m/s Height (m)

0

0.4

Average void fraction (-)

800 ml, 0.032 m/s 1.2 1 0.8 0.6 0.4 0.2 0

0.2

500 ml, 0.100 m/s Height (m)

Height (m) 0

Average void fraction (-)

1.2 1 0.8 0.6 0.4 0.2 0

Height (m)

0.2

1.2 1 0.8 0.6 0.4 0.2 0

Height (m)

0

500 ml, 0.085 m/s

500 ml, 0.056 m/s Height (m)

Height (m)

500 ml, 0.032 m/s 1.2 1 0.8 0.6 0.4 0.2 0

69

1.2 1 0.8 0.6 0.4 0.2 0

1000 ml, 0.093m/s

0

0.2

0.4

0.6

Average void fraction (-)

Ergun Gibilaro GidaspowErgunWenYu GidaspowSchillerNaumann SchillerNaumann SyamlalOBrien WenYu Experiment

Fig. 14. Concentration profiles as solid void fraction, after 20 s of simulation time, 2-D simulations grid size 100 × 1 × 1200.

were beyond the feasible limit, the data points inside the limit confirmed the concentration profile quite well at the top of the bed. Only the experimental points quite close to the feasible limit overpredict the void fraction. The concentration profiles in the simulations were somewhat dependant on the grid size, but the bed height did not show much difference on the used grid. The literature survey shows that the variety of the drag models is large. It also shows that some of the models are derived purely from experiments; some are purely from theory and the rest are a combination of these. The interaction between phases has been formulated either as a drag function (K) or as a drag coefficient (CD). Some of the models have significant differences, which show that the drag models are highly dependant on the system and therefore the CFD modeling cannot be considered as universal as it is normally considered. In 2-dimensional modeling all beds acted as a fluidized bed, where the rise was stable, although in the Ergun model; in the model by Syamal and O'Brien; and in the model by Haider and Levenspiel, when ϕ ≠ 1, the bed did not reach a stable point within the same limit as the other models did. In 3-dimensional modeling these models predicted that the system as a moving bed instead of a fluidized bed. Since the transition regime of this phenomenon is small the eventual result can differ largely from the reality. Acknowledgments The authors like to thank Anssi Salminen for conducting the experiments and Mikko Pihlanko for running part of the CFD work. The authors like to also acknowledge CSC — IT Center for Science Ltd. for the allocation of computational resources and the Finnish Funding Agency for Technology and Innovation for the financial support.

Appendix A. Supplementary data Supplementary data to this article can be found online at doi:10. 1016/j.minpro.2011.12.006. References Asif, M., Ibrahim, A.A., 2002. Minimum fluidization velocity and defluidization behavior of binary-solid liquid-fluidized beds. Powder Technol. 126, 241–254. Atkins, P.W., 1998. Physical chemistry, 6th ed. Oxford University Press, Oxford. Benyahia, S., Syamlal, M., O'Brien, T.J., 2005. Summary of MFIX Equations 2005-4. From URL https://mfix.netl.doe.gov/documentation/MFIXEquations2005-4-4.pdf, 2005. August 2008. Boemer, A., Qi, H., Renz, U., Vasquez, S., Boysan, F., 1995. Eulerian computation of fluidized bed hydrodynamics — a comparison of physical models. ASME FBC Conf. 2, 775–787. Bothe, D., Shirzadi, H., Warnecke, H.-J., 2007. Evaluations of Euler–Euler simulations of bubble columns based on numerical tracer experiments. Chem. Eng. Res. Des. 85 (11), 1491–1496. Boyer, C., Duquenne, A., Wild, G., 2002. Measuring techniques in gas–liquid and gas– liquid-solid reactors. Chem. Eng. Sci. 57, 3185–3215. Burke, S.P., Plummer, W.B., 1928. Gas flow through packed columns. Ind. Eng. Chem. 20, 1169–1200. Carman, P.C., 1937. Fluid flow through granular beds. Trans. Inst. Chem. Eng. (London) 15, 150–166. Dallavalle, J.M., 1948. Micromeritics: the Technology of Fine Particles, 2nd ed. Pitman, London. Drew, D.A., 1983. Mathematical modeling of two-phase flow. Annu. Rev. Fluid Mech. 15 (1), 261–291. Enwald, H., Peirano, E., Almstedt, A.E., 1996. Eulerian two-phase flow theory applied to fluidization. Int. J. Multiphase Flow 22, 21–66 Suppl.. Ergun, S., 1952. Fluid flow through packed columns. Chem. Eng. Prog. 48 (2), 89–94. Garside, J., Al-Dibouni, M.R., 1977. Velocity-Voidage Relationships for Fluidization in Solid-Liquid Systems. Ind. Eng. Chem. Process Des. Dev. 16, 206–214. Gibilaro, L.G., Di Felice, R., Waldram, S.P., 1985. Generalized friction factor and drag coefficient for fluid-particle interactions. Chem. Eng. Sci. 40 (10), 1817–1823. Gidaspow, D., 1994. Multiphase Flow and Fluidization, 1st ed. Academic Press, New York. 467 p.

70

O. Visuri et al. / International Journal of Mineral Processing 104-105 (2012) 58–70

Gidaspow, D., Ettehadieh, B., 1983. Fuidization in two-dimensional beds with a jet. 2. Hydrodynamic modeling. Ind. Eng. Chem. Fundam. 22 (2), 193–201. Goldstein, S., 1929. The steady flow of viscous fluid past a fixed spherical obstacle at small Reynolds numbers. Proc. R. Soc. London, Ser. A 123, 225–235. Gosman, A.D., Issa, R.I., Lekakou, C., Looney, M.K., Politis, S., 1992. Multidimensional modelling of turbulent two-phase flows in stirred vessels. AICHE J. 38 (12), 1946–1956. Haider, A., Levenspiel, O., 1989. Drag coefficient and terminal velocity of spherical and nonspherical particles. Powder Technol. 58, 63–70. Kozeny, J., 1927. Über kapilläre Leitung des Wassers in Boden. Sitzungsber. Akad. Wiss. Wien Math.Naturwiss. Kl. 136, 271–306 Abt.2a. Lundberg, J., Halvorsen, B.M., 2008. A review of some existing drag models describing the interaction between phases in a bubbling fluidized bed. Proc. 49th Scand. Conf. Simulation and Modeling. Oslo University College, Oslo, Norway. 7–8 October. Macdonald, I.F., El-Sayed, M., Mow, K., Dullien, F.A.L., 1979. Flow through porous media — the Ergun equation revisited. Ind. Eng. Chem. Fundam. 18, 199–208. OpenCFD Ltd., 2009. OpenFOAM — The Open Source CFD Toolbox. http://www. openfoam.com/, 2009, accessed 03.02.2011. Oseen, C.W., 1913. Über den Giiltigkeitsbereich der Stokesschen Widerstandsformel. Ark. Mat. Astron. Fys. 9 (16), 1–15. Panneerselvam, R., Savithri, S., Surender, G.D., 2007. CFD based investigations on hydrodynamics and energy dissipation due to solid motion in liquid fluidised bed. Chem. Eng. J. 132, 159–171. Ramesh, P.L.N., Raajenthiren, M., 2010. A Review of Some Existing Drag Models Describing the Interaction between the Solid Gasenous Phases in a CFB. Int. J. Eng. Sci. and Tech. (IJEST) 2 (5), 1047–1051.

Reddy, R.K., Joshi, J.B., 2009. CFD modeling of solid–liquid fluidized beds of mono and binary particle mixtures. Chem. Eng. Sci. 64, 3641–3658. Richardson, J.F., Zaki, W.N., 1954. Sedimentation and fluidization: Part I. Trans. Inst. Chem. Eng. 32, 35–53. Rowe, P.N., 1961. Drag forces in a hydraulic model of a fluidized bed — part II. Trans. Inst. Chem. Eng. 39, 175–180. Rusche, H, 2002. Computational Fluid Dynamics of dispersed two-phase flows at high volume fraction. PhD thesis Imperial College, London, UK. Schiller, L., Naumann, A., 1933. Über die grundlegenden berechnungen bei der schwerkraftauf-bereitung. Z. Ver. Dtsch. Ing. 77 (12), 318–320. Silva, L.F.L.R., Damian, R.B., Lage, P.L.C., 2008. Implementation and analysis of numerical solution of the population balance equation in CFD packages. Comput. Chem. Eng. 32, 2933–2945. Swamee, P.K., Ohja, C.S.P., 1991. Drag coefficient and fall velocity of nonspherical particles. J. Hydraul. Eng. 117, 660–667. Syamlal, M., O’Brien, T.J., 1987. The Derivation of a Drag Coefficient Formula from Velocity-Voidage Correlations. From URL: https://mfix.netl.doe.gov/open_citations/ citation_page.php?citation_id=70, accessed 12.01.2012 Wachem, B. van, 2000. Derivation, Implementation, and Validation of Computer Simulation Models for Gas-Solid Fluidized Beds. PhD thesis Delft University of Technology, The Netherlands. Wen, C.Y., Yu, Y.H., 1966. Mechanics of fluidization: AIChe Symposium Series, 62, pp. 100–111.