**Computational Fluid Dynamics (CFD) and Discrete Element Method (DEM) Applied to Centrifuges**

Xiana Romani Fernandez, Lars Egmont Spelter and Hermann Nirschl *Karlsruhe Institute of Technology, Campus South, Institute for Mechanical Process Engineering and Mechanics, Karlsruhe Germany* 

### **1. Introduction**

96 Applied Computational Fluid Dynamics

[44] H. Kurioka, Y. Oka, H. Satoh, H. Kuwanna, O. Sugawa, Properties of the Plume and

[45] O. Vauquelin, O. Megret, Smoke extraction experiments in case of fire in a tunnel. Fire

[46] Ingason, H., Correlation between temperatures and oxygen measurements in a tunnel

[47] O. Vauquelin, Y. Wu, Influence of tunnel width on longitudinal smoke control. Fire

[48] P. J. Woodburn, R. E. Britter, CFD Simulation of a Tunnel Fire - part two. Fire Safety

[49] G. B. Grant, S. F. Jagger, C. J. Lea, Fires in Tunnels. Philosophical Transactions:

[50] S. R. Lee, H. S. Ryou, A numerical study on smoke movement in longitudinal

[51] R. O. Carvel, A. N. Beard, P. W. Jowitt, The influence of longitudinal ventilation

[52] V. K. R. Kodur, L. A. Bisby, M. F. Green, Experimental evaluation of the fire behaviour

[53] F. Wald, et. al., Experimental behaviour of a steel structure under natural fire. Fire

[54] F. Lestari, et. al., An alternative method for fire smoke toxicity assessment using human

ventilation tunnel fires for different aspect ratio. Building and Environment, 2006.

systems on fires in tunnels. Tunneling and Underground Space Technology, 2001.

of insulated fibre-reinforced-polymer-strengthened reinforced concrete columns.

Methematical, Physical, Engineering Sciences, 1998: p. 2873.

Engineering, 2001(546): p. 151 - 156.

flow. Fire Safety Journal, 2007. 42: p. 75.

Safety Journal, 2002. 37: p. 525.

Safety Journal, 2006. 41: p. 420.

Fire Safety Journal, 2006. 41: p. 547.

lung cells. Fire Safety Journal, 2006. 41: p. 605.

Safety Journal, 2006. 41: p. 509.

Journal, 1996. 26: p. 63.

41: p. 719.

16: p. 3.

Near Fire Source in horizontally long and narrow spaces. Journal of Construction

The simulation of fluid flow in e.g. air channels, pipes, porous media and turbines has been the emphasis of many research projects and optimization processes in the academic and industrial community during the last years. A key advantage of the calculation of fluid flow through a given geometry when compared to the experimental determination is the efficient identification of disadvantages of the current state and the fast cycle of change in geometry and analysis of its effects. This saves time and costs and furthermore opens up new fields of research and development. There are, however, few drawbacks of the new optimisation and research methods. The numerically obtained results need validation before it can be stated whether they are correct or misrepresent the true state. In some cases, such as the optimisation of existing systems, the experience of the user may be sufficient to evaluate the calculated flow pattern. In the case of the prediction of flow in new geometries or process conditions, the numerical result can be validated only by measuring the flow pattern at representative process conditions (flow velocity, pressure, and temperature). For nonmoving geometries, many measurements have been conducted in order to determine the flow behaviour in various geometries and for different process conditions. These measurements are readily available for the evaluation of new simulations. However, only little research has been conducted to evaluate the flow conditions in centrifuges, especially at high rotational speeds. The reason for the lack of experimental data lies in the difficult conditions for the measurement in the centrifugal field with high circumferential velocities and pressure gradients. For low rotational speeds, the flow patterns have been determined by an electrolytic technique (Glinka, 1983) and visualised by adding ink to the flowing liquid (Bass, 1959; 1959; 1962).

The emphasis of the work presented is the evaluation of the efficiency of the acceleration geometry and the prediction of the axial flow profile in the centrifuge. In centrifugation technology, the prediction of the separation efficiency is often restricted to the cut size, the smallest particle size that can settle on the bowl wall at certain operating conditions. In a centrifuge, the highest velocities occur in tangential direction but, due to the throughput, a secondary flow in the axial direction appears. The tangential velocity creates the centrifugal force acting on the particles and the flow in the axial direction determines the residence time

Computational Fluid Dynamics (CFD) and

calculated flow pattern.

**2. Solid bowl centrifuges** 

other apparatus in solid-liquid separation.

Fig. 1. Tubular bowl centrifuge (CEPA Z41G)

Discrete Element Method (DEM) Applied to Centrifuges 99

conditions, which is of high importance for the design and optimization of centrifuges and

This work focuses on the prediction and validation of the results of flow patterns in solid bowl centrifuges of different kind. The centrifuges operate within a range of 1000 up to 40000 rpm and with throughputs between 0.1 and 6 l/min. The article highlights the issues that will occur if flow behaviour has to be calculated in a rotating confinement using commercial software. Detailed results show the variation of the flow patterns achieved using different models and the possibilities to run a sanity check for the

The centrifugal separation of particles in a suspension is one of the most common problems appearing in industrial processes such as waste water, mining and mineral processing, biotechnology, solid-fuel, pharmaceutical, chemical and food industry. Solid bowl centrifuges are widely used in the pharmaceutical and fine chemicals industry. The tubular bowl centrifuge has been providing the basis for all centrifuge designs. A scheme of a

typical state of the art tubular bowl centrifuge is shown in Figure 1.

of particles in the bowl. Thus it is necessary to know the flow patterns in a centrifuge in order to calculate its separation capability. Often the predicted cut size of a centrifuge diverges from the real separation efficiency due to the fact that the assumptions of the analytical models are not strictly valid for industrial process conditions. The deviation is caused by the complex shape and internal assemblies of industrial solid bowl centrifuges, which create complex flow patterns. Hence estimating these flow patterns, as well as measuring them, is complicated.

In rotating machinery, the fluid flow is mostly turbulent, but in centrifuges the degree of turbulence is difficult to determine. Furthermore it may vary throughout the different zones of the geometry. For example, at the inlet of a continuous centrifuge, high shear rates and intensive mixing occur due to the high tangential velocity gradient. Hence a solver that incorporates equations to solve turbulent flow has to be used to obtain accurate predictions. Behind the inlet zone, the liquid flows undisturbed until the continuous phase leaves the centrifuge via an overflow weir or another discharge system. In this zone, the degree of turbulence depends on the geometry of the outlet system and the velocity gradients at the outlet. This example shows how complicated it is to gain the basis for the selection of an appropriate simulation method. Thus it is advisable to choose a model that is as precise as possible, but with an affordable computational demand. Various models have been developed to predict turbulent flow, but the codes were initially not designed to solve the Navier-Stokes equation in a rotating reference frame. Modifications were made over the years to include rotating systems and phenomena, e.g. the modified k-ε model "k-ε-RNG". Since then, only little work has been conducted to evaluate the accuracy of the various turbulence models when applied to centrifuges. The work presented compares several models such as the k-ε RNG, the k-w, the Reynolds Stress Model (RSM) and the large-eddy simulation (LES) for predicting flow pattern in centrifuges. Furthermore the computational demand is evaluated and compared for different mesh sizes and types. The nomographs presented help to estimate the computational time for a given problem and computer system depending on the model and mesh.

Since with increasing fill level the sediment interacts with the flow pattern and reduces the capacity of the centrifuge, the sediment build-up must be included in the calculation of the separation efficiency and the determination of the flow pattern for tubular centrifuges. For calculation of the settling particles, diverse approaches are available. In this study, a Lagrange formulation has been chosen for calculation of the particle trajectories. In order to investigate the sediment build-up, particle–particle interactions must be considered. For this purpose, it is appropriate to apply the Discrete Element Method (DEM) (Cundall & Strack, 1979), originally developed for calculating the flow of bulk granular materials and then extended to the entire particle processing technology. This method accounts for particleparticle and particle-wall interactions by means of various contact models. However this method is limited because it does not account for hydrodynamic forces (drag, lift, torque), which play a vital role in e.g. centrifuges. Hence the simulation of the multiphase flow in centrifuges requires a coupling of the DEM with computational fluid dynamics (CFD). The flow pattern and particle trajectories must be calculated in an alternating scheme. First, the flow pattern is determined using CFD and afterwards, the motion of the particles is calculated using DEM adding the hydrodynamic forces to the contact forces.

The advantage of the numerical simulation of the complex multiphase flow in centrifuges is that the flow, the particle trajectories and their deposition are represented accurately. It allows a detailed description of the flow and particle behaviour for various operating conditions, which is of high importance for the design and optimization of centrifuges and other apparatus in solid-liquid separation.

This work focuses on the prediction and validation of the results of flow patterns in solid bowl centrifuges of different kind. The centrifuges operate within a range of 1000 up to 40000 rpm and with throughputs between 0.1 and 6 l/min. The article highlights the issues that will occur if flow behaviour has to be calculated in a rotating confinement using commercial software. Detailed results show the variation of the flow patterns achieved using different models and the possibilities to run a sanity check for the calculated flow pattern.

### **2. Solid bowl centrifuges**

98 Applied Computational Fluid Dynamics

of particles in the bowl. Thus it is necessary to know the flow patterns in a centrifuge in order to calculate its separation capability. Often the predicted cut size of a centrifuge diverges from the real separation efficiency due to the fact that the assumptions of the analytical models are not strictly valid for industrial process conditions. The deviation is caused by the complex shape and internal assemblies of industrial solid bowl centrifuges, which create complex flow patterns. Hence estimating these flow patterns, as well as

In rotating machinery, the fluid flow is mostly turbulent, but in centrifuges the degree of turbulence is difficult to determine. Furthermore it may vary throughout the different zones of the geometry. For example, at the inlet of a continuous centrifuge, high shear rates and intensive mixing occur due to the high tangential velocity gradient. Hence a solver that incorporates equations to solve turbulent flow has to be used to obtain accurate predictions. Behind the inlet zone, the liquid flows undisturbed until the continuous phase leaves the centrifuge via an overflow weir or another discharge system. In this zone, the degree of turbulence depends on the geometry of the outlet system and the velocity gradients at the outlet. This example shows how complicated it is to gain the basis for the selection of an appropriate simulation method. Thus it is advisable to choose a model that is as precise as possible, but with an affordable computational demand. Various models have been developed to predict turbulent flow, but the codes were initially not designed to solve the Navier-Stokes equation in a rotating reference frame. Modifications were made over the years to include rotating systems and phenomena, e.g. the modified k-ε model "k-ε-RNG". Since then, only little work has been conducted to evaluate the accuracy of the various turbulence models when applied to centrifuges. The work presented compares several models such as the k-ε RNG, the k-w, the Reynolds Stress Model (RSM) and the large-eddy simulation (LES) for predicting flow pattern in centrifuges. Furthermore the computational demand is evaluated and compared for different mesh sizes and types. The nomographs presented help to estimate the computational time for a given problem and computer

Since with increasing fill level the sediment interacts with the flow pattern and reduces the capacity of the centrifuge, the sediment build-up must be included in the calculation of the separation efficiency and the determination of the flow pattern for tubular centrifuges. For calculation of the settling particles, diverse approaches are available. In this study, a Lagrange formulation has been chosen for calculation of the particle trajectories. In order to investigate the sediment build-up, particle–particle interactions must be considered. For this purpose, it is appropriate to apply the Discrete Element Method (DEM) (Cundall & Strack, 1979), originally developed for calculating the flow of bulk granular materials and then extended to the entire particle processing technology. This method accounts for particleparticle and particle-wall interactions by means of various contact models. However this method is limited because it does not account for hydrodynamic forces (drag, lift, torque), which play a vital role in e.g. centrifuges. Hence the simulation of the multiphase flow in centrifuges requires a coupling of the DEM with computational fluid dynamics (CFD). The flow pattern and particle trajectories must be calculated in an alternating scheme. First, the flow pattern is determined using CFD and afterwards, the motion of the particles is

calculated using DEM adding the hydrodynamic forces to the contact forces.

The advantage of the numerical simulation of the complex multiphase flow in centrifuges is that the flow, the particle trajectories and their deposition are represented accurately. It allows a detailed description of the flow and particle behaviour for various operating

measuring them, is complicated.

system depending on the model and mesh.

The centrifugal separation of particles in a suspension is one of the most common problems appearing in industrial processes such as waste water, mining and mineral processing, biotechnology, solid-fuel, pharmaceutical, chemical and food industry. Solid bowl centrifuges are widely used in the pharmaceutical and fine chemicals industry. The tubular bowl centrifuge has been providing the basis for all centrifuge designs. A scheme of a typical state of the art tubular bowl centrifuge is shown in Figure 1.

Fig. 1. Tubular bowl centrifuge (CEPA Z41G)

Computational Fluid Dynamics (CFD) and

pattern for tubular centrifuges.

Fig. 2. Flow pattern in solid bowl centrifuges

**2.1 Flow pattern** 

Discrete Element Method (DEM) Applied to Centrifuges 101

and a rotor without any settled particles that may interact with the flow pattern. In fact, the length L is reduced by up to 30 % due to the effects of the inlet. With increasing fill level, the sediment interacts with the flow pattern and reduces the free length L so that the cut size of the centrifuge increases (Spelter et al., 2010). Therefore the sediment build-up has to be included in the calculation of the separation efficiency and the determination of the flow

There are two limiting cases for the flow pattern in solid bowl centrifuges. Either the entire volume is passed by the feed or the main axial fluid flow takes place in a thin layer on top of a stagnant but rotating pool. The two different flow patterns are schematically shown in Figure 2. When a boundary layer flow is present, the geometric constant rbowl in Equation 1

When starting building centrifuges and evaluating the separation performance, it was assumed that the proportion of the volume in which significant axial fluid velocities occur, further referred to as active volume Vac or area Aac, is nearly as high as the entire volume of the centrifuge, Vac=100 % (Horanyi & Nemeth, 1971). Other researchers realised that the active volume of the centrifuge is significantly smaller than the capacity of the centrifuge (Vac<<100 %) and developed the boundary layer theory (Bass, 1959). Today, most researchers who work with solid bowl centrifuges support the boundary layer theory.

Recent studies by the authors have confirmed the validity of this theory for industrial solid bowl and tubular bowl centrifuges (Romani Fernandez & Nirschl, 2009; Spelter & Nirschl, 2010; Spelter et al., 2010). Nevertheless it is important to pay attention to the acceleration geometry of the centrifuge and the influence of the settled particles on the flow pattern. In many centrifuges, the feed enters the rotor via a cone or a jet hitting a distributor plate. The

has to be replaced by the radius where the stagnant pool begins (rboundary layer).

Starting with a rotating cylinder, modifications were made, leading to the disc stack centrifuge, screen scroll centrifuge and various other types. The basic design still has its applications. Due to the lean cylinder, high rotational speeds are possible. One of the first high speed centrifuges was the Sharples superfuge (Saunders, 1948; Taylor, 1946) that has been modified and adapted by various companies. Applications range from laboratory purposes to the production of pharmaceuticals in an industrial scale. New centrifuge designs were developed allowing the extension of the techniques to the industrial manufacture of new products such as vaccines (Anderson et al., 1969; Brantley et al., 1970; Perardi & Anderson, 1970; Perardi et al., 1969). The application of most high-speed centrifuges is limited by the solids capacity of the rotor. Due to the lean design, the solids capacity usually does not exceed a few kilograms. For most products that are processed in such centrifuges, the small capacity does not increase the production costs significantly because the value of the final product exceeds the cost of labour and machinery by an order of magnitude. This is the reason why most research has been conducted in the field of process optimization, e.g. in sterilisation, product recovery, and operator safety. Optimisation of the flow pattern in the tubular centrifuge has not yet been the focus of research.

The advantages of this type of centrifuge are the high centrifugal numbers and the simple geometry as compared to e.g. disc stack centrifuges. Tubular bowl centrifuges reach centrifugal numbers up to 40000 g. Their stability is due to their relatively high length to diameter aspect ratio in the range of approximately 4 to 7. There are other types of solid bowl basket centrifuges with conical or tubular bowls but smaller length to diameter aspect ratios, usually below 0.75 (Leung, 1998) such as overflow centrifuges. Sometimes radial or vertical compartments as well as blades are included to improve the separation efficiency.

Up to now, the flow pattern in existing centrifuges has been analysed and compared to some theories describing the fluid flow in rotating geometries (Glinka, 1983; Golovko, 1969; Gösele, 1968; 1974; Schaeflinger & Stibi, 1991; Schaeflinger, 1991; Trawinski, 1959). Using the determined axial flow pattern, the authors developed various analytical models for calculation of the cut size. Since it was not possible to measure or predict the tangential velocity, corrections were made to take into account insufficient tangential acceleration. These empirical corrections were often only valid for a specific centrifuge.

Moreover, the cut size depends not only on the tangential and axial flow but also on the density difference Δρ between the solids and the liquid, the viscosity of the suspension μ and the geometry of the rotor as stated in Equation (1). The geometric parameters are the length of the rotor L, the radius of the overflow weir rweir and the radius of the bowl rbowl:

$$\chi = \sqrt{\frac{18 \cdot \mu \cdot \text{v}\_{\text{ax}} \cdot \ln\left(\frac{\mathbf{r}\_{\text{bowl}}}{\mathbf{r}\_{\text{weir}}}\right)}{\Delta \rho \cdot \text{L}}} \cdot \frac{2 \cdot \pi \cdot \text{n}}{60} \tag{1}$$

Modifications were made because the observed separation efficiency of the test rig often was not in good agreement with simple theoretical approaches. These modifications included certain levels of turbulence, boundary layer flow and inhomogeneous flow in the acceleration and discharge zones (Bass, 1959; 1959; 1962; Sokolow, 1971; Zubkov & Golovko, 1968; 1969). The different flow patterns in centrifuges are explained in detail in Chapter 2.1. The assumption that the entire length L is available for the deposition of particles is only valid for an ideal tangential acceleration of the feed so that no turbulent inlet area is present and a rotor without any settled particles that may interact with the flow pattern. In fact, the length L is reduced by up to 30 % due to the effects of the inlet. With increasing fill level, the sediment interacts with the flow pattern and reduces the free length L so that the cut size of the centrifuge increases (Spelter et al., 2010). Therefore the sediment build-up has to be included in the calculation of the separation efficiency and the determination of the flow pattern for tubular centrifuges.

### **2.1 Flow pattern**

100 Applied Computational Fluid Dynamics

Starting with a rotating cylinder, modifications were made, leading to the disc stack centrifuge, screen scroll centrifuge and various other types. The basic design still has its applications. Due to the lean cylinder, high rotational speeds are possible. One of the first high speed centrifuges was the Sharples superfuge (Saunders, 1948; Taylor, 1946) that has been modified and adapted by various companies. Applications range from laboratory purposes to the production of pharmaceuticals in an industrial scale. New centrifuge designs were developed allowing the extension of the techniques to the industrial manufacture of new products such as vaccines (Anderson et al., 1969; Brantley et al., 1970; Perardi & Anderson, 1970; Perardi et al., 1969). The application of most high-speed centrifuges is limited by the solids capacity of the rotor. Due to the lean design, the solids capacity usually does not exceed a few kilograms. For most products that are processed in such centrifuges, the small capacity does not increase the production costs significantly because the value of the final product exceeds the cost of labour and machinery by an order of magnitude. This is the reason why most research has been conducted in the field of process optimization, e.g. in sterilisation, product recovery, and operator safety. Optimisation of the flow pattern in the

The advantages of this type of centrifuge are the high centrifugal numbers and the simple geometry as compared to e.g. disc stack centrifuges. Tubular bowl centrifuges reach centrifugal numbers up to 40000 g. Their stability is due to their relatively high length to diameter aspect ratio in the range of approximately 4 to 7. There are other types of solid bowl basket centrifuges with conical or tubular bowls but smaller length to diameter aspect ratios, usually below 0.75 (Leung, 1998) such as overflow centrifuges. Sometimes radial or vertical compartments as well as blades are included to improve the separation efficiency. Up to now, the flow pattern in existing centrifuges has been analysed and compared to some theories describing the fluid flow in rotating geometries (Glinka, 1983; Golovko, 1969; Gösele, 1968; 1974; Schaeflinger & Stibi, 1991; Schaeflinger, 1991; Trawinski, 1959). Using the determined axial flow pattern, the authors developed various analytical models for calculation of the cut size. Since it was not possible to measure or predict the tangential velocity, corrections were made to take into account insufficient tangential acceleration.

Moreover, the cut size depends not only on the tangential and axial flow but also on the density difference Δρ between the solids and the liquid, the viscosity of the suspension μ and the geometry of the rotor as stated in Equation (1). The geometric parameters are the length of the rotor L, the radius of the overflow weir rweir and the radius of the bowl rbowl:

> L r

Modifications were made because the observed separation efficiency of the test rig often was not in good agreement with simple theoretical approaches. These modifications included certain levels of turbulence, boundary layer flow and inhomogeneous flow in the acceleration and discharge zones (Bass, 1959; 1959; 1962; Sokolow, 1971; Zubkov & Golovko, 1968; 1969). The different flow patterns in centrifuges are explained in detail in Chapter 2.1. The assumption that the entire length L is available for the deposition of particles is only valid for an ideal tangential acceleration of the feed so that no turbulent inlet area is present

 

weir bowl

 

<sup>r</sup> <sup>18</sup> <sup>v</sup> ln

ax

60 2 n

(1)

These empirical corrections were often only valid for a specific centrifuge.

x

tubular centrifuge has not yet been the focus of research.

There are two limiting cases for the flow pattern in solid bowl centrifuges. Either the entire volume is passed by the feed or the main axial fluid flow takes place in a thin layer on top of a stagnant but rotating pool. The two different flow patterns are schematically shown in Figure 2. When a boundary layer flow is present, the geometric constant rbowl in Equation 1 has to be replaced by the radius where the stagnant pool begins (rboundary layer).

When starting building centrifuges and evaluating the separation performance, it was assumed that the proportion of the volume in which significant axial fluid velocities occur, further referred to as active volume Vac or area Aac, is nearly as high as the entire volume of the centrifuge, Vac=100 % (Horanyi & Nemeth, 1971). Other researchers realised that the active volume of the centrifuge is significantly smaller than the capacity of the centrifuge (Vac<<100 %) and developed the boundary layer theory (Bass, 1959). Today, most researchers who work with solid bowl centrifuges support the boundary layer theory.

Fig. 2. Flow pattern in solid bowl centrifuges

Recent studies by the authors have confirmed the validity of this theory for industrial solid bowl and tubular bowl centrifuges (Romani Fernandez & Nirschl, 2009; Spelter & Nirschl, 2010; Spelter et al., 2010). Nevertheless it is important to pay attention to the acceleration geometry of the centrifuge and the influence of the settled particles on the flow pattern. In many centrifuges, the feed enters the rotor via a cone or a jet hitting a distributor plate. The

Computational Fluid Dynamics (CFD) and

cylinder is 200 mm.

without hindrance.

setup; b) outlet for tandem setup

**3. Mathematical formulation of the fluid** 

mathematical methodology used to calculate these parameters.

Discrete Element Method (DEM) Applied to Centrifuges 103

inner and outer cylinder at a radius of 42 mm at the top. The inlet is marked with an arrow and the height of the inlet is further referred to as H=0 mm. Both cylinders and the feed pipes rotate with the same revolutions per minute. There is no air in the system as the void is entirely filled with water. The liquid leaves the centrifuge at two outlet bores at a radius of 36 mm. In the overflow setup, only the outer cylinder is used and thus an air core is present. The liquid leaves the centrifuge via the two outlet bores at 36 mm and via two additional outlet bores at a radius of 30 mm, which adjust the minimum pool level. The two bores at 36 mm can be sealed. The height of the outlet is at H=170 mm, the overall length of the

Depending on the pressure drop at the outlet bores, the water level rises above the radius of the weir. The superelevation is the equilibrium of pressure drop, rotational speed and throughput. With higher rotational speed and lower throughput the elevation is close to zero. With increasing throughput and lower rotational speeds, the centrifuge can be entirely filled. In the test rig, the maximum elevation is limited by a bore (radius 12.5 mm) that is necessary for the measurement of the flow pattern with Laser Doppler Anemometry (Spelter, et al. 2011). Once this radius is reached by the liquid, it overflows the circular weir

Fig. 4. Simplified scheme of the setups of the centrifuge test rig - outer rotor made of carbon fibre reinforced plastic, inner one made of tempered glass - a) additional outlet for overflow

To describe the flow pattern, it is necessary to determine the velocity components in all of the spatial directions considered as well as the pressure and, in the case of compressible flow, the density at any time and position in the domain. This chapter summarizes the

The conservation equations of the fluid mechanics, conservation of mass, momentum and energy, can be obtained by means of a balance around an infinitesimally small fluid element. The continuity equation, which is derived from the mass conservation, describes

liquid is not accelerated to the tangential velocity of the pool surface and thus hits the surface with a certain velocity gradient. The rotating pool accelerates the incoming liquid to the angular velocity. The pressure gradient in the rotating liquid stabilises the flow pattern and reduces the mixing of the feed and the rotating pool in radial direction. Thus a stratification of the liquid in radial layers occurs from the inlet and the liquid flows axially at the liquid-gas interface in a boundary layer. The flow is eased by the reduced friction at the interface between liquid and gas. The flow pattern will be different if the liquid is not fed to the pool surface or if intensive mixing occurs at the inlet. This is the case for some tubular centrifuges, as long as no sediment is present. Figure 3 a) shows the inlet geometry of a tubular bowl centrifuge as it was used in previous studies. The metal structure accelerates the incoming liquid to a certain rotational speed and significant mixing of incoming and rotating liquid occurs. The rigid body rotation is reached after a few centimetres behind the inlet. Figure 3 b) shows the acceleration geometry after 5 minutes of centrifugation. The settled yeast cells block a significant section of the inlet geometry and thus impede intensive mixing and acceleration. This has been observed for various minerals and biological products (Stahl et al., 2008) and is detectable by measuring the residence time and separation efficiency in the centrifuge at different fill levels. The growth of the sediment will change the flow pattern significantly if it reaches the active area. This may not be the case for industrial scale solid bowl centrifuges but has to be included in the calculations of the separation efficiency in tubular bowl centrifuges.

Fig. 3. Inlet geometry of tubular bowl centrifuge – a) empty rotor; b) with sediment after 5 min of centrifugation at 15000 rpm; c) as "b" but at 40000 rpm

#### **2.2 Rotor geometry**

Other inlet systems have been analysed in a recent research project (Spelter, et al. 2011). The liquid is accelerated by rotating inlet pipes and the outlet of these pipes is below the pool surface. This minimises the lag of the tangential velocity of the incoming liquid significantly and the trapping of air in the fluid. This feeding system may be used in existing centrifuges in order to enhance the separation efficiency by reducing the axial section where the tangential velocity gradient reduces the effective centrifugal force on the settling particles. Figure 4 shows a simplified scheme of the rotor of the test rig. The test rig is operated with two different rotor setups. Either the centrifuge runs using both cylinders, further referred to as tandem setup, or without the inner glass cylinder, further referred to as overflow setup. The tandem setup is similar to the one applied in high-performance centrifuges used for the production of vaccines. In these centrifuges, the outer rotor is often made of titanium and the inner one of plastic, e.g. Noryl®. The liquid enters the cylinder halfway between

liquid is not accelerated to the tangential velocity of the pool surface and thus hits the surface with a certain velocity gradient. The rotating pool accelerates the incoming liquid to the angular velocity. The pressure gradient in the rotating liquid stabilises the flow pattern and reduces the mixing of the feed and the rotating pool in radial direction. Thus a stratification of the liquid in radial layers occurs from the inlet and the liquid flows axially at the liquid-gas interface in a boundary layer. The flow is eased by the reduced friction at the interface between liquid and gas. The flow pattern will be different if the liquid is not fed to the pool surface or if intensive mixing occurs at the inlet. This is the case for some tubular centrifuges, as long as no sediment is present. Figure 3 a) shows the inlet geometry of a tubular bowl centrifuge as it was used in previous studies. The metal structure accelerates the incoming liquid to a certain rotational speed and significant mixing of incoming and rotating liquid occurs. The rigid body rotation is reached after a few centimetres behind the inlet. Figure 3 b) shows the acceleration geometry after 5 minutes of centrifugation. The settled yeast cells block a significant section of the inlet geometry and thus impede intensive mixing and acceleration. This has been observed for various minerals and biological products (Stahl et al., 2008) and is detectable by measuring the residence time and separation efficiency in the centrifuge at different fill levels. The growth of the sediment will change the flow pattern significantly if it reaches the active area. This may not be the case for industrial scale solid bowl centrifuges but has to be included in the calculations of the

Fig. 3. Inlet geometry of tubular bowl centrifuge – a) empty rotor; b) with sediment after 5

Other inlet systems have been analysed in a recent research project (Spelter, et al. 2011). The liquid is accelerated by rotating inlet pipes and the outlet of these pipes is below the pool surface. This minimises the lag of the tangential velocity of the incoming liquid significantly and the trapping of air in the fluid. This feeding system may be used in existing centrifuges in order to enhance the separation efficiency by reducing the axial section where the tangential velocity gradient reduces the effective centrifugal force on the settling particles. Figure 4 shows a simplified scheme of the rotor of the test rig. The test rig is operated with two different rotor setups. Either the centrifuge runs using both cylinders, further referred to as tandem setup, or without the inner glass cylinder, further referred to as overflow setup. The tandem setup is similar to the one applied in high-performance centrifuges used for the production of vaccines. In these centrifuges, the outer rotor is often made of titanium and the inner one of plastic, e.g. Noryl®. The liquid enters the cylinder halfway between

separation efficiency in tubular bowl centrifuges.

min of centrifugation at 15000 rpm; c) as "b" but at 40000 rpm

**2.2 Rotor geometry** 

inner and outer cylinder at a radius of 42 mm at the top. The inlet is marked with an arrow and the height of the inlet is further referred to as H=0 mm. Both cylinders and the feed pipes rotate with the same revolutions per minute. There is no air in the system as the void is entirely filled with water. The liquid leaves the centrifuge at two outlet bores at a radius of 36 mm. In the overflow setup, only the outer cylinder is used and thus an air core is present. The liquid leaves the centrifuge via the two outlet bores at 36 mm and via two additional outlet bores at a radius of 30 mm, which adjust the minimum pool level. The two bores at 36 mm can be sealed. The height of the outlet is at H=170 mm, the overall length of the cylinder is 200 mm.

Depending on the pressure drop at the outlet bores, the water level rises above the radius of the weir. The superelevation is the equilibrium of pressure drop, rotational speed and throughput. With higher rotational speed and lower throughput the elevation is close to zero. With increasing throughput and lower rotational speeds, the centrifuge can be entirely filled. In the test rig, the maximum elevation is limited by a bore (radius 12.5 mm) that is necessary for the measurement of the flow pattern with Laser Doppler Anemometry (Spelter, et al. 2011). Once this radius is reached by the liquid, it overflows the circular weir without hindrance.

Fig. 4. Simplified scheme of the setups of the centrifuge test rig - outer rotor made of carbon fibre reinforced plastic, inner one made of tempered glass - a) additional outlet for overflow setup; b) outlet for tandem setup

#### **3. Mathematical formulation of the fluid**

To describe the flow pattern, it is necessary to determine the velocity components in all of the spatial directions considered as well as the pressure and, in the case of compressible flow, the density at any time and position in the domain. This chapter summarizes the mathematical methodology used to calculate these parameters.

The conservation equations of the fluid mechanics, conservation of mass, momentum and energy, can be obtained by means of a balance around an infinitesimally small fluid element. The continuity equation, which is derived from the mass conservation, describes

Computational Fluid Dynamics (CFD) and

**3.2 Turbulence models** 

fluctuation:

Reynolds stresses tensor τt:

Discrete Element Method (DEM) Applied to Centrifuges 105

Using the VOF method, a single momentum equation, the Reynolds-averaged Navier-Stokes equation, is solved throughout the domain, and the resulting velocity field is shared among the phases. The dependency on the volume fraction is implemented by using volume-

v v p v F

<sup>2</sup>

 <sup>n</sup>

 <sup>n</sup>

Due to high velocities and complex flow areas, most of the flow patterns of interest in process engineering, aviation, and automotive engineering become unstable above certain Reynolds numbers and, thus turbulences emerge. Turbulent eddies appear at very different time and length scales and are always three-dimensional. The kinetic energy stored in the rotational movement of the large eddies is passed to smaller eddies, until they disappear by energy dissipation. A turbulent flow is characterized by a chaotic and random fluctuation of the flow properties velocity and pressure. The conservation equations are able to capture the physics of the fluctuation motion, but in order to solve all the vortex structures of a turbulent flow, an extremely fine grid and very small time steps would be necessary. The memory and the computing power required for the socalled Direct Numerical Simulation (DNS) can only be provided by supercomputers. For that reason, turbulence is usually modelled with different mathematical approaches. A large-eddy simulation (LES) performs a filtering of the fluctuating flow quantities. The small eddies have a common behaviour and are nearly isotropic, while large eddies are more anisotropic and their behaviour is strongly influenced by the geometry of the domain, the boundary conditions and the body forces. Thus larger eddies are directly calculated, while fine structures are modelled. As a result, the computational demand decreases regarding to the DNS but it is nonetheless high in comparison with the

The flow variables can be represented as the sum of a time-averaged size and an additional

The substitution of these expressions in the momentum conservation equation for incompressible flows and the subsequent temporal averaging yield to the Reynoldsaveraged conservation equations. In the modified continuity equation, the velocity component is just replaced by the time-averaged variable. In the Reynolds-averaged Navier-Stokes equation (RANS), the fluctuation velocities appear in an additional tensor, the

u u u' ; v v v' ; w w w' ; p p p' . (9)

t

<sup>q</sup> <sup>1</sup> <sup>q</sup> <sup>q</sup> (7)

<sup>q</sup> <sup>1</sup> <sup>q</sup> <sup>q</sup> (8)

(6)

averaged values for the density and viscosity , as explained in Equations 6 - 8.

t v

approaches modelling all the turbulence scales.

 

the temporal change of mass in the volume element as the net rate of flow mass into the element across its boundaries. In the case of the incompressible flow considered here, the continuity equation simplifies to

$$
\nabla \cdot \vec{\mathbf{v}} = 0 \; . \tag{2}
$$

In an analogous way, the momentum equation can be determined based on the conservation of momentum in all of the spatial directions considered. Newton's second law states that the rate of change of momentum of a fluid element equals the sum of the forces acting on the volume element. These forces can be surface forces, such as viscous forces and shear and normal stresses, or mass forces, such as gravity and centrifugal forces. For incompressible and Newtonian fluids in turbulent regime, the resulting momentum conservation equation follows to be

$$\rho \left( \frac{\partial \vec{\mathbf{v}}}{\partial \mathbf{t}} + (\vec{\mathbf{v}} \cdot \nabla) \vec{\mathbf{v}} \right) = -\nabla \mathbf{p} + \mu \nabla^2 \vec{\mathbf{v}} + \nabla \cdot \mathbf{r}\_{\mathbf{t}} + \vec{\mathbf{F}} \tag{3}$$

where τt represents the tensor of turbulences and F other additional forces that can be included expressed as volumetric forces. All terms from Equation (3), the Reynoldsaveraged Navier-Stokes equation (RANS), are discretized and calculated for each volume cell with the exception of τt, which is modelled. The RANS equation together with the continuity equation form an equation system of four differential equations with four unknown variables, the flow velocities u, v, w in x, y and z direction, respectively and the pressure p, which can be numerically calculated. This equation will be used for simulation of the centrifugal flow when only one phase is considered. For the cases where the liquid pool rotates around an air core, a multiphase approach must be used in order to simulate both, the air and the liquid. The presence of solid particles is ignored for flow simulation purposes, which is an acceptable assumption for low solid concentrations of the suspensions, as is the case in the present study.

#### **3.1 Two-phase approach**

The Volume of Fluid method (VOF) (Hirt & Nichols, 1981) is an interface tracking method that simulates the gas-liquid multiphase flow. This method is a simple and efficient formulation designed to track the interface of two phases that do not interpenetrate with a relatively small number of interfaces. Hence it is often used to simulate the gas-liquid multiphase flow in industrial devices (Brennan, 2003; Brennan et al., 2007; Li et al., 1999; Mousavian & Najafi, 2009). VOF introduces a variable which takes values from zero to one and represents the volume fraction of one of the phases in each cell. A continuity conservation equation is solved for each phase

$$\frac{\partial}{\partial t} \left( \alpha\_{\mathbf{q}} \rho\_{\mathbf{q}} \right) + \nabla \cdot \left( \alpha\_{\mathbf{q}} \rho\_{\mathbf{q}} \vec{\mathbf{v}}\_{\mathbf{q}} \right) = 0 \tag{4}$$

Furthermore the volume fraction of each phase q in any cell has to obey

$$\sum\_{\mathbf{q}=1}^{n} \mathbf{a}\_{\mathbf{q}} = \mathbf{1} \,. \tag{5}$$

Using the VOF method, a single momentum equation, the Reynolds-averaged Navier-Stokes equation, is solved throughout the domain, and the resulting velocity field is shared among the phases. The dependency on the volume fraction is implemented by using volumeaveraged values for the density and viscosity , as explained in Equations 6 - 8.

$$\overline{\rho} \left( \frac{\partial \vec{\mathbf{v}}}{\partial t} + (\vec{\mathbf{v}} \cdot \nabla) \vec{\mathbf{v}} \right) = -\nabla \mathbf{p} + \overline{\mu} \, \nabla^2 \vec{\mathbf{v}} + \nabla \cdot \mathbf{\tau}\_t + \vec{\mathbf{F}} \tag{6}$$

$$\overline{\rho} = \sum\_{\mathbf{q}=1}^{n} \alpha\_{\mathbf{q}} \ \rho\_{\mathbf{q}} \tag{7}$$

$$
\overline{\mu} = \sum\_{\mathbf{q}=1}^{n} \alpha\_{\mathbf{q}} \; \mu\_{\mathbf{q}} \tag{8}
$$

#### **3.2 Turbulence models**

104 Applied Computational Fluid Dynamics

the temporal change of mass in the volume element as the net rate of flow mass into the element across its boundaries. In the case of the incompressible flow considered here, the

In an analogous way, the momentum equation can be determined based on the conservation of momentum in all of the spatial directions considered. Newton's second law states that the rate of change of momentum of a fluid element equals the sum of the forces acting on the volume element. These forces can be surface forces, such as viscous forces and shear and normal stresses, or mass forces, such as gravity and centrifugal forces. For incompressible and Newtonian fluids in turbulent regime, the resulting momentum conservation equation

v v p v F

<sup>2</sup>

where τt represents the tensor of turbulences and F other additional forces that can be included expressed as volumetric forces. All terms from Equation (3), the Reynoldsaveraged Navier-Stokes equation (RANS), are discretized and calculated for each volume cell with the exception of τt, which is modelled. The RANS equation together with the continuity equation form an equation system of four differential equations with four unknown variables, the flow velocities u, v, w in x, y and z direction, respectively and the pressure p, which can be numerically calculated. This equation will be used for simulation of the centrifugal flow when only one phase is considered. For the cases where the liquid pool rotates around an air core, a multiphase approach must be used in order to simulate both, the air and the liquid. The presence of solid particles is ignored for flow simulation purposes, which is an acceptable assumption for low solid concentrations of the

The Volume of Fluid method (VOF) (Hirt & Nichols, 1981) is an interface tracking method that simulates the gas-liquid multiphase flow. This method is a simple and efficient formulation designed to track the interface of two phases that do not interpenetrate with a relatively small number of interfaces. Hence it is often used to simulate the gas-liquid multiphase flow in industrial devices (Brennan, 2003; Brennan et al., 2007; Li et al., 1999; Mousavian & Najafi, 2009). VOF introduces a variable which takes values from zero to one and represents the volume fraction of one of the phases in each cell. A continuity

v 0

<sup>1</sup> <sup>n</sup> q 1 <sup>q</sup>

.

. (5)

q<sup>q</sup> q<sup>q</sup> <sup>q</sup> 

t

Furthermore the volume fraction of each phase q in any cell has to obey

t v

suspensions, as is the case in the present study.

conservation equation is solved for each phase

**3.1 Two-phase approach** 

 

<sup>v</sup> <sup>0</sup> . (2)

t

,

(3)

(4)

continuity equation simplifies to

follows to be

Due to high velocities and complex flow areas, most of the flow patterns of interest in process engineering, aviation, and automotive engineering become unstable above certain Reynolds numbers and, thus turbulences emerge. Turbulent eddies appear at very different time and length scales and are always three-dimensional. The kinetic energy stored in the rotational movement of the large eddies is passed to smaller eddies, until they disappear by energy dissipation. A turbulent flow is characterized by a chaotic and random fluctuation of the flow properties velocity and pressure. The conservation equations are able to capture the physics of the fluctuation motion, but in order to solve all the vortex structures of a turbulent flow, an extremely fine grid and very small time steps would be necessary. The memory and the computing power required for the socalled Direct Numerical Simulation (DNS) can only be provided by supercomputers. For that reason, turbulence is usually modelled with different mathematical approaches. A large-eddy simulation (LES) performs a filtering of the fluctuating flow quantities. The small eddies have a common behaviour and are nearly isotropic, while large eddies are more anisotropic and their behaviour is strongly influenced by the geometry of the domain, the boundary conditions and the body forces. Thus larger eddies are directly calculated, while fine structures are modelled. As a result, the computational demand decreases regarding to the DNS but it is nonetheless high in comparison with the approaches modelling all the turbulence scales.

The flow variables can be represented as the sum of a time-averaged size and an additional fluctuation:

$$\mathbf{u} = \overline{\mathbf{u}} + \mathbf{u}'; \qquad \mathbf{v} = \overline{\mathbf{v}} + \mathbf{v}'; \qquad \mathbf{w} = \overline{\mathbf{w}} + \mathbf{w}'; \qquad \mathbf{p} = \overline{\mathbf{p}} + \mathbf{p}'.\tag{9}$$

The substitution of these expressions in the momentum conservation equation for incompressible flows and the subsequent temporal averaging yield to the Reynoldsaveraged conservation equations. In the modified continuity equation, the velocity component is just replaced by the time-averaged variable. In the Reynolds-averaged Navier-Stokes equation (RANS), the fluctuation velocities appear in an additional tensor, the Reynolds stresses tensor τt:

Computational Fluid Dynamics (CFD) and

& Nirschl, 2009; Spelter & Nirschl, 2010).

were used in this work.

with the density and viscosity of water.

Discrete Element Method (DEM) Applied to Centrifuges 107

standard k-ε model. The k-ε RNG model has been succeeding for previous simulations of the flow in centrifuges by other groups and in our own research group (Romani Fernandez

The k-ω model (Wilcox, 1988), like the k-ε model, uses two equations to describe turbulence. Instead of the dissipation rate ε, the specific dissipation or turbulent frequency ω=k/ε is introduced in this model. The turbulent viscosity is then given by Equation (13). The accuracy of the model depends on the undisturbed velocity of the fluid outside the boundary layer which is subject to large fluctuations. This dependence can lead to significant errors in the calculation (Bardina et al., 1997). Due to the strong dissipation of the turbulence on the wall, no special treatment of the boundary layer at the wall is necessary.

> 

Comparing the turbulence models, the k-ω model, as a two-equation model, has a similar range of strengths and weaknesses as the k- ε model. RSM is complex, but it is generally accepted as the simplest kind of model with the potential to describe all mean flow properties and Reynolds stresses without case-by-case adjustment. RSM is by no means as comprehensively validated as the k-ε model and due to the high computational cost it is not as widely used in industrial flow simulations as the k-ε model. LES, due to the inherent unsteady nature, is much more computationally expensive than the k-ε and k-ω models. However, compared to RSM it was proved that it requires only twice the computational demand of RSM (Ferziger, 1977) as cited in (Versteeg & Malalasekera, 2007). This is not a big difference taking into account the solution accuracy and the ability of LES to resolve certain time-dependent features. In order to investigate the applicability of the different turbulence models to the simulation of the flow pattern in centrifugal field, different approaches, the k-ε RNG model, the k-ω model, the RSM, and the LES,

In order to simplify the simulation, the air is considered as an incompressible gas. This is a reasonable assumption for the operation conditions of atmospheric pressure and for a nontemperature dependent problem. The fluid is defined as an incompressible Newtonian one

The inlet velocity is set to match the desired volume flux. Using a turbulence model, also the turbulence quantities of the inlet flow must be specified. The input of fixed values for k and ε or k and ω are difficult to estimate. Thus these parameters are usually obtained as a function of the turbulence intensity and a characteristic length which must be set at the boundaries (Paschedag, 2004). For the multiphase simulations, a value of one was given to the volume fraction of water at the inlet. A static pressure of 101325 Pa (1 atmosphere) is set as ambient condition for the outlet, defined as a pressure outlet. Sometimes, there is a backflow through the pressure outlet in order to obey the mass conservation. Therefore it is recommended to set realistic conditions for a possible backflow to avoid convergence

It is possible to reduce the computational demand by dividing the geometry into periodic sections. The flow pattern is calculated in one of these segments with periodic conditions at

<sup>k</sup> t (13)

This feature is desired for the exact calculation of the flow near the wall.

**3.3 Boundary conditions, discretization schemata and solver** 

problems. For the cases presented, only backflow of air is allowed.

$$\mathbf{r}\_t = -\rho \begin{pmatrix} \overline{\mathbf{u'} \mathbf{u'}} & \overline{\mathbf{u'} \mathbf{v'}} & \overline{\mathbf{u'} \mathbf{w'}}\\ \overline{\mathbf{v'} \mathbf{u'}} & \overline{\mathbf{w'} \mathbf{v'}} & \overline{\mathbf{v'} \mathbf{w'}}\\ \overline{\mathbf{w'} \mathbf{u'}} & \overline{\mathbf{w'} \mathbf{v'}} & \overline{\mathbf{w'} \mathbf{w'}} \end{pmatrix}. \tag{10}$$

The description of a turbulent flow with the RANS equation requires modelling of these tensor terms. The most complex classical turbulence model is the Reynolds Stress Model, RSM, (Launder et al., 1975), which uses transport equations to model each of the elements of the stress tensor and in addition an equation for the dissipation rate of the turbulences. This means that five additional transport equations are required in two-dimensional flows and seven in three-dimensional flow. This way, the directional effects of the Reynolds stress field can be taken into account. This model is recommended for flows with complex strain fields, such as highly swirling flows, or significant body forces. However, the computational cost of RSM is often unaffordable.

An effective method to reduce the modelling effort is to apply the Boussinesq approach, Equation (11). This approach introduces the eddy viscosity or turbulent viscosity μt, which represents the momentum and energy transport by diffusion by the eddies or turbulent fluctuations. Thus the Reynolds stresses can be calculated as the viscous stresses and only this new variable μt must be modelled.

$$-\rho \,\overline{\mathbf{u'}\,\mathbf{v'}} = -\rho \,\overline{\mathbf{v'}\,\mathbf{u'}} = \mu\_t \left(\frac{\partial \overline{\mathbf{u}}}{\partial \mathbf{y}} + \frac{\partial \overline{\mathbf{v}}}{\partial \mathbf{x}}\right) \tag{11}$$

One of the most common turbulence models for the determination of μt is the k-ε model (Launder & Spalding, 1974). This model introduces two variables to calculate μt, the turbulent kinetic energy k, and its dissipation rate ε:

$$
\mu\_{\rm t} = \rho \mathbf{C}\_{\mu} \frac{\mathbf{k}^2}{\mathbf{c}},
\tag{12}
$$

where Cμ represents a semi-empirical dimensionless constant. The temporal and spatial changes of both variables are described with a transport equation which includes not just the convection and diffusion but also the creation of turbulence and its dissipation by using several source and sink terms. This advantage justifies the additional expense of the transport equation turbulence models compared to algebraic turbulence models, which presume that the turbulence only depends on local conditions. The disadvantage of the k-ε model is that it assumes the turbulences to be isotropic, which is not strictly true in most of the flows. This model, well established and most widely validated, has already been used to simulate the flow in cyclones and centrifuges; although in its development this model assumes a fully turbulent flow, which is only partially applicable to centrifuges.

The k-ε renormalization group model, k-ε RNG, (Yakhot & Orszag, 1986) is an extension of the standard k-ε model that takes into account the effect of swirl on the turbulence by means of an extra source term in the transport equations. Thus the k-ε RNG model exhibits a higher accuracy for swirling flows. The RNG procedure systematically removes the small scales of motion from the governing equations by expressing their effects in terms of larger scale motions and a modified viscosity (Versteeg & Malalasekera, 2007). While the standard k- ε model is appropriate for high-Reynolds-number flows, the RNG theory accounts for low-Reynolds-number effects. Thus this model is more reliable for a wider class of flows than the

w' u' w' v' w' w' v' u' v' v' v' w' u' u' u' v' u' w'

The description of a turbulent flow with the RANS equation requires modelling of these tensor terms. The most complex classical turbulence model is the Reynolds Stress Model, RSM, (Launder et al., 1975), which uses transport equations to model each of the elements of the stress tensor and in addition an equation for the dissipation rate of the turbulences. This means that five additional transport equations are required in two-dimensional flows and seven in three-dimensional flow. This way, the directional effects of the Reynolds stress field can be taken into account. This model is recommended for flows with complex strain fields, such as highly swirling flows, or significant body forces. However, the computational cost of

An effective method to reduce the modelling effort is to apply the Boussinesq approach, Equation (11). This approach introduces the eddy viscosity or turbulent viscosity μt, which represents the momentum and energy transport by diffusion by the eddies or turbulent fluctuations. Thus the Reynolds stresses can be calculated as the viscous stresses and only

One of the most common turbulence models for the determination of μt is the k-ε model (Launder & Spalding, 1974). This model introduces two variables to calculate μt, the

where Cμ represents a semi-empirical dimensionless constant. The temporal and spatial changes of both variables are described with a transport equation which includes not just the convection and diffusion but also the creation of turbulence and its dissipation by using several source and sink terms. This advantage justifies the additional expense of the transport equation turbulence models compared to algebraic turbulence models, which presume that the turbulence only depends on local conditions. The disadvantage of the k-ε model is that it assumes the turbulences to be isotropic, which is not strictly true in most of the flows. This model, well established and most widely validated, has already been used to simulate the flow in cyclones and centrifuges; although in its development this model

The k-ε renormalization group model, k-ε RNG, (Yakhot & Orszag, 1986) is an extension of the standard k-ε model that takes into account the effect of swirl on the turbulence by means of an extra source term in the transport equations. Thus the k-ε RNG model exhibits a higher accuracy for swirling flows. The RNG procedure systematically removes the small scales of motion from the governing equations by expressing their effects in terms of larger scale motions and a modified viscosity (Versteeg & Malalasekera, 2007). While the standard k- ε model is appropriate for high-Reynolds-number flows, the RNG theory accounts for low-Reynolds-number effects. Thus this model is more reliable for a wider class of flows than the

t

assumes a fully turbulent flow, which is only partially applicable to centrifuges.

 

RSM is often unaffordable.

this new variable μt must be modelled.

turbulent kinetic energy k, and its dissipation rate ε:

 

<sup>t</sup> . (10)

 

u' v' v' u' t (11)

<sup>k</sup> <sup>C</sup> , (12)

 

x v

 

2

y u standard k-ε model. The k-ε RNG model has been succeeding for previous simulations of the flow in centrifuges by other groups and in our own research group (Romani Fernandez & Nirschl, 2009; Spelter & Nirschl, 2010).

The k-ω model (Wilcox, 1988), like the k-ε model, uses two equations to describe turbulence. Instead of the dissipation rate ε, the specific dissipation or turbulent frequency ω=k/ε is introduced in this model. The turbulent viscosity is then given by Equation (13). The accuracy of the model depends on the undisturbed velocity of the fluid outside the boundary layer which is subject to large fluctuations. This dependence can lead to significant errors in the calculation (Bardina et al., 1997). Due to the strong dissipation of the turbulence on the wall, no special treatment of the boundary layer at the wall is necessary. This feature is desired for the exact calculation of the flow near the wall.

$$
\mu\_t = \rho \frac{\mathbf{k}}{\varpi} \tag{13}
$$

Comparing the turbulence models, the k-ω model, as a two-equation model, has a similar range of strengths and weaknesses as the k- ε model. RSM is complex, but it is generally accepted as the simplest kind of model with the potential to describe all mean flow properties and Reynolds stresses without case-by-case adjustment. RSM is by no means as comprehensively validated as the k-ε model and due to the high computational cost it is not as widely used in industrial flow simulations as the k-ε model. LES, due to the inherent unsteady nature, is much more computationally expensive than the k-ε and k-ω models. However, compared to RSM it was proved that it requires only twice the computational demand of RSM (Ferziger, 1977) as cited in (Versteeg & Malalasekera, 2007). This is not a big difference taking into account the solution accuracy and the ability of LES to resolve certain time-dependent features. In order to investigate the applicability of the different turbulence models to the simulation of the flow pattern in centrifugal field, different approaches, the k-ε RNG model, the k-ω model, the RSM, and the LES, were used in this work.

#### **3.3 Boundary conditions, discretization schemata and solver**

In order to simplify the simulation, the air is considered as an incompressible gas. This is a reasonable assumption for the operation conditions of atmospheric pressure and for a nontemperature dependent problem. The fluid is defined as an incompressible Newtonian one with the density and viscosity of water.

The inlet velocity is set to match the desired volume flux. Using a turbulence model, also the turbulence quantities of the inlet flow must be specified. The input of fixed values for k and ε or k and ω are difficult to estimate. Thus these parameters are usually obtained as a function of the turbulence intensity and a characteristic length which must be set at the boundaries (Paschedag, 2004). For the multiphase simulations, a value of one was given to the volume fraction of water at the inlet. A static pressure of 101325 Pa (1 atmosphere) is set as ambient condition for the outlet, defined as a pressure outlet. Sometimes, there is a backflow through the pressure outlet in order to obey the mass conservation. Therefore it is recommended to set realistic conditions for a possible backflow to avoid convergence problems. For the cases presented, only backflow of air is allowed.

It is possible to reduce the computational demand by dividing the geometry into periodic sections. The flow pattern is calculated in one of these segments with periodic conditions at

Computational Fluid Dynamics (CFD) and

condition, for the pressure correction.

**4. Simulated flow in solid bowl centrifuges** 

technique for the determination of flow velocities.

number of 0.25 was applied for the volume fraction.

Discrete Element Method (DEM) Applied to Centrifuges 109

fraction is refined with respect to the general time step based on the maximal Courant number allowed for this variable. The Courant number, also called Courant-Froudel-Lewy number (CFL), Equation 15, is a dimensionless number that compares the time step of the calculation to the characteristic time needed by a fluid element to cross a cell. In some explicit time discretization methods, the criterion CFL<1 (Oertel & Laurien, 2003) is imposed to achieve stability in the numerical calculation. An explicit schema limited by a Courant

xcell /vi

(15)

<sup>t</sup> CFL 

The simulations were performed with the commercial software ANSYS FLUENT with the versions 6.3.26, 12.0.16, 12.1.2 and 13.0, which are comparable. An overview of the type of solver used by ANSYS can be found in (Kelecy, 2008). The solver used is an unsteady, segregated pressure-based solver. The segregated pressure-based algorithm solves the conservation equations sequentially. Because these equations are coupled, the solution loop must be carried out iteratively in order to obtain a converged numerical solution. Pressure and velocity were coupled either with the SIMPLE algorithm, Semi-Implicit Method for Pressure Linked Equations, or with the PISO schema, which stands for Pressure Implicit with Splitting of Operators and is available for unsteady calculations. The PISO pressurevelocity coupling algorithm has one predictor step and two corrector steps, extending the pressure correction procedure of the commonly used SIMPLE algorithm with a further corrector step (Versteeg & Malalasekera, 2007). The SIMPLE algorithm is a guess-andcorrect procedure which uses a Poisson's equation, developed from the continuity

The previously published simulations of a tubular bowl centrifuge were carried out twodimensionally (Spelter & Nirschl, 2010); hence it was not possible to include the real acceleration geometry in the model. For the work presented, all simulations were run threedimensionally, including the acceleration and discharge geometry, and were compared with experimental data. The flow pattern was predicted for a tubular bowl centrifuge with a feeding system similar to the one used in a recent published work (Spelter, et al. 2011). The experimental data has been obtained by Laser Doppler Anemometry, a non-invasive

**4.1 Prediction of the flow pattern (k-**e **RNG model, transient solution) - tandem setup**  Figure 5 shows the tangential velocity profiles for four different rotational speeds at a constant throughput of 1 l/min. There is no significant deviation between the predicted tangential velocity and the rigid body rotation (RBR). The measurement of the tangential velocity for various rotational speeds and throughputs at different heights was in good agreement with the rigid body rotation as well (Spelter, et al. 2011). Thus it can be stated that the simulation is in good agreement with the experimentally determined tangential velocity profiles. Therefore it is possible to optimise existing acceleration geometries by using a feeding system similar to the one proposed. By calculating the present state and modifying the geometry, the changes in the acceleration efficiency can be evaluated by

the intersecting planes. The periodic surfaces were defined with a periodic boundary condition which uses the flow conditions at one of the periodic surfaces as its cells were the adjacent cells of the other periodic plane. The rotational periodicity must be taken into account and the rotation axis must be defined. Due to the rotational nature of the system, it is helpful to use a rotating reference frame to perform the calculations. The side and top walls of the bowl, as well as the inlet accelerator, all defined with no-slip condition, rotate with the same angular velocity as the reference frame.

Another problem emerging from the use of turbulence models is the treatment of turbulent quantities at the walls. Despite the turbulent flow, there always exists a laminar sublayer at the wall followed by a transient region until the flow can be considered as turbulent. A common approach is to create a first layer of cells next to the wall with a width including the entire viscous sublayer and the transition layer. In this layer, standard wall functions (Paschedag, 2004), most widely used for industrial flows, can be applied. The wall functions define values for the velocity field and the turbulent quantities suitable as boundary conditions for the solution in the second cells layer, where the turbulence is fully developed. This approach is particularly suitable for flows with larger spatial extent in which the thickness of the boundary layer at the wall is small compared to the extension of the whole domain. To choose the thickness of the first layer, the criterion of the y+ value (Equation 14) is often used. The value of y+ is a function of the first layer thickness y, the shear stress at the wall τw, and the physical properties of the fluid density and viscosity μ. Values of y+ between 30 and 100 allow the use of standard wall functions. For details concerning the exact formulation of the y+ and standard wall functions, the reader is referred to (ANSYS, 2009) which is based on (Launder & Spalding, 1974).

$$\mathbf{y}^+ = \frac{\rho \text{ y}}{\mu} \sqrt{\frac{|\mathbf{r}\_\mathbf{w}|}{\rho}} \tag{14}$$

The k-ε model and the RSM need these wall functions in order to calculate the turbulent properties near the wall. The k-ω model does not require wall-damping functions, because the value of turbulence kinetic energy at the wall is set to zero and for the turbulent frequency ω a hyperbolic variation at the near-wall grid point is applied (Wilcox, 1988).

The discretization schema used for the velocity and the turbulent variables is set to firstorder upwind, due to the convective nature of the flow in centrifugal field and to avoid convergence problems. For the volume fraction, the Geo-Reconstruct method (ANSYS, 2006) is applied in order to obtain a sharp interface between both phases. For the interface between two fluids, this method assumes a linear slope, which is calculated with a piecewise linear approach. The variable pressure was discretized with the PRESTO! (PREssure STaggering Option) scheme recommended for high speed rotating flows (ANSYS, 2009). If velocities and pressure values are both defined at the centres of the cells, a highly non-uniform pressure field can act like a uniform field in the discretized momentum equations (Versteeg & Malalasekera, 2007). A solution to this problem is to use a staggered grid to calculate the pressure via a discrete continuity balance. In the staggered grid, the values of pressure in the centre are known and these are the values at the interfaces in the normal grid.

The temporal discretization of the transport properties was performed using a first-order implicit method. However, using the VOF model, the time step for the variable volume

the intersecting planes. The periodic surfaces were defined with a periodic boundary condition which uses the flow conditions at one of the periodic surfaces as its cells were the adjacent cells of the other periodic plane. The rotational periodicity must be taken into account and the rotation axis must be defined. Due to the rotational nature of the system, it is helpful to use a rotating reference frame to perform the calculations. The side and top walls of the bowl, as well as the inlet accelerator, all defined with no-slip condition, rotate

Another problem emerging from the use of turbulence models is the treatment of turbulent quantities at the walls. Despite the turbulent flow, there always exists a laminar sublayer at the wall followed by a transient region until the flow can be considered as turbulent. A common approach is to create a first layer of cells next to the wall with a width including the entire viscous sublayer and the transition layer. In this layer, standard wall functions (Paschedag, 2004), most widely used for industrial flows, can be applied. The wall functions define values for the velocity field and the turbulent quantities suitable as boundary conditions for the solution in the second cells layer, where the turbulence is fully developed. This approach is particularly suitable for flows with larger spatial extent in which the thickness of the boundary layer at the wall is small compared to the extension of the whole domain. To choose the thickness of the first layer, the criterion of the y+ value (Equation 14) is often used. The value of y+ is a function of the first layer thickness y, the shear stress at the wall τw, and the physical properties of the fluid density and viscosity μ. Values of y+ between 30 and 100 allow the use of standard wall functions. For details concerning the exact formulation of the y+ and standard wall functions, the reader is referred to (ANSYS, 2009) which is based on (Launder & Spalding,

> 

<sup>y</sup> <sup>w</sup> <sup>y</sup> (14)

The k-ε model and the RSM need these wall functions in order to calculate the turbulent properties near the wall. The k-ω model does not require wall-damping functions, because the value of turbulence kinetic energy at the wall is set to zero and for the turbulent frequency ω a hyperbolic variation at the near-wall grid point is applied (Wilcox, 1988). The discretization schema used for the velocity and the turbulent variables is set to firstorder upwind, due to the convective nature of the flow in centrifugal field and to avoid convergence problems. For the volume fraction, the Geo-Reconstruct method (ANSYS, 2006) is applied in order to obtain a sharp interface between both phases. For the interface between two fluids, this method assumes a linear slope, which is calculated with a piecewise linear approach. The variable pressure was discretized with the PRESTO! (PREssure STaggering Option) scheme recommended for high speed rotating flows (ANSYS, 2009). If velocities and pressure values are both defined at the centres of the cells, a highly non-uniform pressure field can act like a uniform field in the discretized momentum equations (Versteeg & Malalasekera, 2007). A solution to this problem is to use a staggered grid to calculate the pressure via a discrete continuity balance. In the staggered grid, the values of pressure in the centre are known and these are the values at

The temporal discretization of the transport properties was performed using a first-order implicit method. However, using the VOF model, the time step for the variable volume

with the same angular velocity as the reference frame.

1974).

the interfaces in the normal grid.

fraction is refined with respect to the general time step based on the maximal Courant number allowed for this variable. The Courant number, also called Courant-Froudel-Lewy number (CFL), Equation 15, is a dimensionless number that compares the time step of the calculation to the characteristic time needed by a fluid element to cross a cell. In some explicit time discretization methods, the criterion CFL<1 (Oertel & Laurien, 2003) is imposed to achieve stability in the numerical calculation. An explicit schema limited by a Courant number of 0.25 was applied for the volume fraction.

$$\text{CFL} = \frac{\text{At}}{\Delta \mathbf{x}\_{\text{cell}} / \mathbf{v}\_{\text{i}}} \tag{15}$$

The simulations were performed with the commercial software ANSYS FLUENT with the versions 6.3.26, 12.0.16, 12.1.2 and 13.0, which are comparable. An overview of the type of solver used by ANSYS can be found in (Kelecy, 2008). The solver used is an unsteady, segregated pressure-based solver. The segregated pressure-based algorithm solves the conservation equations sequentially. Because these equations are coupled, the solution loop must be carried out iteratively in order to obtain a converged numerical solution. Pressure and velocity were coupled either with the SIMPLE algorithm, Semi-Implicit Method for Pressure Linked Equations, or with the PISO schema, which stands for Pressure Implicit with Splitting of Operators and is available for unsteady calculations. The PISO pressurevelocity coupling algorithm has one predictor step and two corrector steps, extending the pressure correction procedure of the commonly used SIMPLE algorithm with a further corrector step (Versteeg & Malalasekera, 2007). The SIMPLE algorithm is a guess-andcorrect procedure which uses a Poisson's equation, developed from the continuity condition, for the pressure correction.

### **4. Simulated flow in solid bowl centrifuges**

The previously published simulations of a tubular bowl centrifuge were carried out twodimensionally (Spelter & Nirschl, 2010); hence it was not possible to include the real acceleration geometry in the model. For the work presented, all simulations were run threedimensionally, including the acceleration and discharge geometry, and were compared with experimental data. The flow pattern was predicted for a tubular bowl centrifuge with a feeding system similar to the one used in a recent published work (Spelter, et al. 2011). The experimental data has been obtained by Laser Doppler Anemometry, a non-invasive technique for the determination of flow velocities.

#### **4.1 Prediction of the flow pattern (k-**e **RNG model, transient solution) - tandem setup**

Figure 5 shows the tangential velocity profiles for four different rotational speeds at a constant throughput of 1 l/min. There is no significant deviation between the predicted tangential velocity and the rigid body rotation (RBR). The measurement of the tangential velocity for various rotational speeds and throughputs at different heights was in good agreement with the rigid body rotation as well (Spelter, et al. 2011). Thus it can be stated that the simulation is in good agreement with the experimentally determined tangential velocity profiles. Therefore it is possible to optimise existing acceleration geometries by using a feeding system similar to the one proposed. By calculating the present state and modifying the geometry, the changes in the acceleration efficiency can be evaluated by

Computational Fluid Dynamics (CFD) and

comparison to analytical values

Laser Doppler Anemometry as well.

bowl wall.

Discrete Element Method (DEM) Applied to Centrifuges 111

The gauge pressure must be taken into account for dimensioning the rotor. Furthermore it is possible to increase the density of the liquid to model the arbitrary load of a sediment on the

Fig. 6. Pressure distribution, all values in bar: a) 9000 rpm, pressure b) 1000-14000 rpm and

The axial velocity profiles are displayed as contour plots in Figure 7 for 1000 rpm (a) and 3000 rpm (b) at a constant throughput of 1 l/min. The inlet is located at the bottom while the liquid leaves the centrifuge at the top. The incoming liquid enters the centrifuge as a jet that does not widen up significantly towards the outlet. At 1000 rpm, the main axial flow takes place at a small circular cross-section that is not displaced in angular direction on the path towards the outlet bores. With increasing rotational speed, a small angular displacement of the axial flow path is predicted, as shown in Figure 7 b). This effect has been observed by

Fig. 7. Axial velocity in m/s, throughput 1 l/min: a) 1000 rpm b) 3000 rpm

means of CFD. Furthermore the calculation will help in the case of a scale-up of the centrifuge. The acceleration efficiency and the shear in the feed zone can be calculated before a pilot machine or a full scale centrifuge is built. With the knowledge of the flow behaviour, the geometry can be optimised, leading to an improved apparatus with less cost in development and testing.

Fig. 5. Comparison of calculated tangential velocity and rigid body rotation, 1 l/min

The tangential velocity of the fluid creates a pressure in the cylinder in addition to the stress in the rotor created by the rotating cylinder itself. For the case that the rotor is designed as an ideal cylinder, the pressure distribution and thus the tensile stress can be calculated analytically (Sokolow, 1971; Stahl, 2004). However in many industrial centrifuges, the rotor has a non-cylindrical shape. For these designs, the stress in the rotor is determined via finite element methods (FEM), but the arbitrary load by the rotating liquid has to be calculated separately. The geometry of the rotor can be imported from various CAD programs such as ProEngineer, Autodesk Inventor, Catia, etc.. The pressure on the different planes in the bowl is then accessible by calculating the pressure distribution for the highest rotational speed that occurs in the centrifuge. The combination of FEM and CFD leads to a complete virtual development and saves time and costs. The simulations will never substitute completely experimental evaluation of the separation efficiency and mechanical integrity of the centrifuge, but the scale-up factor may be significantly increased. The predicted pressure distributions in the cylinder are shown in Figure 6 for 9000 rpm (a) and for a range of rotational speeds in (b). The outlet is located at the top end of the cylinder. The predicted pressures in Figure 6 b) are compared with the analytically calculated pressure distribution for the different rotational speeds. All values are given in bar.

means of CFD. Furthermore the calculation will help in the case of a scale-up of the centrifuge. The acceleration efficiency and the shear in the feed zone can be calculated before a pilot machine or a full scale centrifuge is built. With the knowledge of the flow behaviour, the geometry can be optimised, leading to an improved apparatus with less cost

Fig. 5. Comparison of calculated tangential velocity and rigid body rotation, 1 l/min

The tangential velocity of the fluid creates a pressure in the cylinder in addition to the stress in the rotor created by the rotating cylinder itself. For the case that the rotor is designed as an ideal cylinder, the pressure distribution and thus the tensile stress can be calculated analytically (Sokolow, 1971; Stahl, 2004). However in many industrial centrifuges, the rotor has a non-cylindrical shape. For these designs, the stress in the rotor is determined via finite element methods (FEM), but the arbitrary load by the rotating liquid has to be calculated separately. The geometry of the rotor can be imported from various CAD programs such as ProEngineer, Autodesk Inventor, Catia, etc.. The pressure on the different planes in the bowl is then accessible by calculating the pressure distribution for the highest rotational speed that occurs in the centrifuge. The combination of FEM and CFD leads to a complete virtual development and saves time and costs. The simulations will never substitute completely experimental evaluation of the separation efficiency and mechanical integrity of the centrifuge, but the scale-up factor may be significantly increased. The predicted pressure distributions in the cylinder are shown in Figure 6 for 9000 rpm (a) and for a range of rotational speeds in (b). The outlet is located at the top end of the cylinder. The predicted pressures in Figure 6 b) are compared with the analytically calculated pressure distribution for the different rotational speeds. All

in development and testing.

values are given in bar.

The gauge pressure must be taken into account for dimensioning the rotor. Furthermore it is possible to increase the density of the liquid to model the arbitrary load of a sediment on the bowl wall.

Fig. 6. Pressure distribution, all values in bar: a) 9000 rpm, pressure b) 1000-14000 rpm and comparison to analytical values

The axial velocity profiles are displayed as contour plots in Figure 7 for 1000 rpm (a) and 3000 rpm (b) at a constant throughput of 1 l/min. The inlet is located at the bottom while the liquid leaves the centrifuge at the top. The incoming liquid enters the centrifuge as a jet that does not widen up significantly towards the outlet. At 1000 rpm, the main axial flow takes place at a small circular cross-section that is not displaced in angular direction on the path towards the outlet bores. With increasing rotational speed, a small angular displacement of the axial flow path is predicted, as shown in Figure 7 b). This effect has been observed by Laser Doppler Anemometry as well.

Fig. 7. Axial velocity in m/s, throughput 1 l/min: a) 1000 rpm b) 3000 rpm

Computational Fluid Dynamics (CFD) and

**4.2 Comparison of different turbulence models** 

inlet and outlet at the angular position of the inlet.

w, RSM and LES turbulence model

3000 rpm

Discrete Element Method (DEM) Applied to Centrifuges 113

Fig. 9. Axial velocity, radius of inlet for different distances from inlet: left: 1000 rpm; right:

It is presumed that the deviations between observed flow pattern and predicted flow behaviour are caused by the turbulence model used. Some authors already compared different models that simplify turbulent effects by correlation factors and additional equations like the k-ε RNG model. They stated that the k-ε RNG model is a reasonable compromise between computational demand and validity of the results (Wardle et al., 2006). Other authors relied on the more extensive Reynolds Stress Model (RSM) or the large-eddy simulation (LES) (Mainza et al., 2006; Wang & Yu, 2006). Nowakowski et al. reviewed the milestones in the application of computational fluid dynamics in hydrocyclones (Nowakowski et al., 2004). The most common models used are the k-ε and recently the RSM and LES model. For the evaluation of the differences in the prediction of the flow behaviour in centrifuges between several turbulence models, the flow pattern was calculated for various rotational speeds and throughputs, using the k-e RNG, k-w, Reynolds Stress Model and large-eddy simulation. A representative comparison of the obtained results by using the different models is shown in Figure 10. The plots illustrate the differences in the calculated tangential and axial velocity halfway between

Fig. 10. Tangential (left) and axial velocity profiles (right), throughput 1 l/min, k-e RNG, k-

However, the predicted axial flow pattern does not depict the experimentally determined flow behaviour accurately. It has been observed that the inlet jet widens within the first 30 % of the rotor (Spelter, et al. 2011) until a parabolic profile is developed that is distributed over the whole circumference. The predicted axial flow profile close to the inlet is in good agreement with the experimental data but the change of flow pattern in direction of the angle of rotation is not correct. The axial flow profiles averaged over the length of the centrifuge between 1000 and 14 000 rpm for two angles of rotation and a constant throughput of 1 l/min are shown in Figure 8. The standard deviation is plotted in both diagrams. The scatter of the data is very high for the angle where the inlet is located, which is 0°. The diagram on the right side of Figure 8 shows the axial flow profiles for an angle of rotation of ±45° displacement to the inlet. There is no significant flow towards the exit of the rotor, only recirculating currents are predicted. The standard deviation of the data is small when compared to the angle of 0°. The origin of the high standard deviation of the data displayed in the diagram on the left in Figure 8 is caused by the change in flow pattern from the feed towards the liquid discharge. The radius of the inlet is higher than the radius of the outlet. Hence the main flow is shifted towards the rotational axis. These effects are displayed in Figure 9 for 1000 rpm (left) and 3000 rpm (right).

Fig. 8. Axial velocity, mean values from throughout the rotor: left: radius of inlet; right: angular position ± 45°

The maximum of the axial velocity is close behind the inlet (H=0mm) at a radius of 42 mm, which corresponds to the inlet radius. Towards the liquid discharge, the maximum shifts to a radius of 37 mm, where the outlet pipes are located. The change in flow pattern is indicated by the arrows in Figure 9. This behaviour was detected for all rotational speeds and different throughputs. The changing flow profile explains the high standard deviation of the averaged profiles shown in the diagram on the left in Figure 8.

The deviation between simulation and experimental data with respect to the flow profile over the full angle of rotation may be explained by the turbulent model used. The modelling does not calculate the turbulences in detail, but estimates the degree of turbulence and thus extenuates the fluctuations at the inlet and outlet. Furthermore the flow between inlet and outlet may not be entirely turbulent, as mentioned in the introduction.

However, the predicted axial flow pattern does not depict the experimentally determined flow behaviour accurately. It has been observed that the inlet jet widens within the first 30 % of the rotor (Spelter, et al. 2011) until a parabolic profile is developed that is distributed over the whole circumference. The predicted axial flow profile close to the inlet is in good agreement with the experimental data but the change of flow pattern in direction of the angle of rotation is not correct. The axial flow profiles averaged over the length of the centrifuge between 1000 and 14 000 rpm for two angles of rotation and a constant throughput of 1 l/min are shown in Figure 8. The standard deviation is plotted in both diagrams. The scatter of the data is very high for the angle where the inlet is located, which is 0°. The diagram on the right side of Figure 8 shows the axial flow profiles for an angle of rotation of ±45° displacement to the inlet. There is no significant flow towards the exit of the rotor, only recirculating currents are predicted. The standard deviation of the data is small when compared to the angle of 0°. The origin of the high standard deviation of the data displayed in the diagram on the left in Figure 8 is caused by the change in flow pattern from the feed towards the liquid discharge. The radius of the inlet is higher than the radius of the outlet. Hence the main flow is shifted towards the rotational axis. These effects are displayed in Figure 9 for 1000 rpm (left) and 3000 rpm

Fig. 8. Axial velocity, mean values from throughout the rotor: left: radius of inlet; right:

of the averaged profiles shown in the diagram on the left in Figure 8.

outlet may not be entirely turbulent, as mentioned in the introduction.

The maximum of the axial velocity is close behind the inlet (H=0mm) at a radius of 42 mm, which corresponds to the inlet radius. Towards the liquid discharge, the maximum shifts to a radius of 37 mm, where the outlet pipes are located. The change in flow pattern is indicated by the arrows in Figure 9. This behaviour was detected for all rotational speeds and different throughputs. The changing flow profile explains the high standard deviation

The deviation between simulation and experimental data with respect to the flow profile over the full angle of rotation may be explained by the turbulent model used. The modelling does not calculate the turbulences in detail, but estimates the degree of turbulence and thus extenuates the fluctuations at the inlet and outlet. Furthermore the flow between inlet and

(right).

angular position ± 45°

Fig. 9. Axial velocity, radius of inlet for different distances from inlet: left: 1000 rpm; right: 3000 rpm

#### **4.2 Comparison of different turbulence models**

It is presumed that the deviations between observed flow pattern and predicted flow behaviour are caused by the turbulence model used. Some authors already compared different models that simplify turbulent effects by correlation factors and additional equations like the k-ε RNG model. They stated that the k-ε RNG model is a reasonable compromise between computational demand and validity of the results (Wardle et al., 2006). Other authors relied on the more extensive Reynolds Stress Model (RSM) or the large-eddy simulation (LES) (Mainza et al., 2006; Wang & Yu, 2006). Nowakowski et al. reviewed the milestones in the application of computational fluid dynamics in hydrocyclones (Nowakowski et al., 2004). The most common models used are the k-ε and recently the RSM and LES model. For the evaluation of the differences in the prediction of the flow behaviour in centrifuges between several turbulence models, the flow pattern was calculated for various rotational speeds and throughputs, using the k-e RNG, k-w, Reynolds Stress Model and large-eddy simulation. A representative comparison of the obtained results by using the different models is shown in Figure 10. The plots illustrate the differences in the calculated tangential and axial velocity halfway between inlet and outlet at the angular position of the inlet.

Fig. 10. Tangential (left) and axial velocity profiles (right), throughput 1 l/min, k-e RNG, kw, RSM and LES turbulence model

Computational Fluid Dynamics (CFD) and

displaced few degrees in the angle of rotation.

observed in the experiments, is formed.

**overflow setup** 

Discrete Element Method (DEM) Applied to Centrifuges 115

at various rotational speeds. The water is fed through an annulus in the two dimensional model and the tangential acceleration takes place via the friction at the wall at both sides of the annulus. According to the prediction, this feeding geometry is sufficient for an acceleration of the incoming liquid up to 9000 rpm. The differences between the cases are distinguishable in the axial flow profiles which are shown in the diagram on the right side of Figure 11. The profiles presented are the values halfway between inlet and outlet. Due to the symmetry condition, the flow is homogeneously distributed along the whole circumference for a specific radius. Therefore the maximum axial velocity is lower than in the three dimensional case due to the conservation of mass. There are recirculating currents next to the inlet in the two dimensional cases which increase in magnitude with rotational speed. In the three dimensional case, recirculating flows are located next to the inlet as well, but

Nevertheless the predicted axial velocity profiles in the two dimensional cases are comparable to the ones obtained in the three dimensional ones in respect to the radial distribution of the flow. The inlet jet remains until the water leaves the centrifuge via the circular weir, so that no plug-shaped profile - resulting in a high active volume -, as

**4.4 Calculation of the flow pattern (k-**e **RNG model, VOF method, transient solution) -** 

The overflow setup was simulated by using the previously described Volume of Fluid method. The centrifuge is filled with water at the beginning of the simulation. As an initial condition, the water is already spinning with the rigid body rotation. This assumption is justifiable because in the experiments, the centrifuge was slowly filled with water and ran for at least 30 minutes in steady state prior to the measurements. The initial condition reduces the necessary time until the steady state is reached significantly. Figure 12 shows

Fig. 12. Profiles of water-air interface for a throughput of 1 l/min: a) 1000 rpm; b) 9000 rpm the distribution of the two phases for 1000 (a) und 9000 rpm (b) for a throughput of 1 l/min. Due to the gravitational and centrifugal force, a conical shape of the water surface is developed at 1000 rpm. The centrifugal number outbalances the gravitational force with

The comparison reveals no significant deviation of the predicted flow pattern. The differences in the tangential velocity are barely detectable. All calculations are in good agreement with the rigid body rotation and thus with the experimentally determined profiles (Spelter, et al. 2011). Using the LES model, the obtained axial flow profile differs from the other models by predicting an earlier shift of the main axial flow towards the radius of the outlet at 36 mm. Because of the high computational demand and instability of the RSM, the simulation with this turbulence model was started with a solution obtained with the k-e RNG model. Nevertheless the differences are small and all calculated axial flow pattern are in disagreement with the measured flow profile. The mixing and expanding of the inlet jet is significantly understated. By not questioning and validating the simulations, the user of the CFD software is led to false conclusions which may result in low performance of the centrifuge due to disadvantageous design. The classical evaluation of the simulation does not indicate a computational error. The mass balance between the inlet and the outlet is fulfilled as stated in Table 1, the residuals are sufficiently small (10-3 to 10-5) and the flow patterns do not change with increasing flow time. Furthermore the predictions are reproducible with different turbulence models. All these facts suggest an accurate calculation. Only the comparison of the main velocities with experimentally determined data shows that the prediction misrepresents the true case in respect to the axial flow profile.


Table 1. Deviation of mass balance dV for different turbulence models in % at 3000 rpm and 1 l/min

#### **4.3 Comparison of two and three dimensional modelling**

The analysis of the acceleration geometry is only possible in three dimensional modelling, because, as mentioned in the introduction, the inlet geometry cannot be correctly reproduced in the two dimensional case. The results presented in chapter 4.1 and 4.2 show a short-circuit flow between inlet and outlet. The comparison between the results obtained when calculating three and two dimensional are shown in Figure 11. The predicted tangential velocity profiles are in good agreement with the rigid body rotation for both cases

Fig. 11. Tangential and axial velocity profiles, throughput 1 l/min, k-e RNG, k-w, RSM and LES turbulence model

The comparison reveals no significant deviation of the predicted flow pattern. The differences in the tangential velocity are barely detectable. All calculations are in good agreement with the rigid body rotation and thus with the experimentally determined profiles (Spelter, et al. 2011). Using the LES model, the obtained axial flow profile differs from the other models by predicting an earlier shift of the main axial flow towards the radius of the outlet at 36 mm. Because of the high computational demand and instability of the RSM, the simulation with this turbulence model was started with a solution obtained with the k-e RNG model. Nevertheless the differences are small and all calculated axial flow pattern are in disagreement with the measured flow profile. The mixing and expanding of the inlet jet is significantly understated. By not questioning and validating the simulations, the user of the CFD software is led to false conclusions which may result in low performance of the centrifuge due to disadvantageous design. The classical evaluation of the simulation does not indicate a computational error. The mass balance between the inlet and the outlet is fulfilled as stated in Table 1, the residuals are sufficiently small (10-3 to 10-5) and the flow patterns do not change with increasing flow time. Furthermore the predictions are reproducible with different turbulence models. All these facts suggest an accurate calculation. Only the comparison of the main velocities with experimentally determined data shows that the prediction misrepresents the true case in

**Model** k-e RNG k-w RSM LES **% dV** 0.0395 0.0015 0.0022 0.0172 Table 1. Deviation of mass balance dV for different turbulence models in % at 3000 rpm and

The analysis of the acceleration geometry is only possible in three dimensional modelling, because, as mentioned in the introduction, the inlet geometry cannot be correctly reproduced in the two dimensional case. The results presented in chapter 4.1 and 4.2 show a short-circuit flow between inlet and outlet. The comparison between the results obtained when calculating three and two dimensional are shown in Figure 11. The predicted tangential velocity profiles are in good agreement with the rigid body rotation for both cases

Fig. 11. Tangential and axial velocity profiles, throughput 1 l/min, k-e RNG, k-w, RSM and

**4.3 Comparison of two and three dimensional modelling** 

respect to the axial flow profile.

1 l/min

LES turbulence model

at various rotational speeds. The water is fed through an annulus in the two dimensional model and the tangential acceleration takes place via the friction at the wall at both sides of the annulus. According to the prediction, this feeding geometry is sufficient for an acceleration of the incoming liquid up to 9000 rpm. The differences between the cases are distinguishable in the axial flow profiles which are shown in the diagram on the right side of Figure 11. The profiles presented are the values halfway between inlet and outlet. Due to the symmetry condition, the flow is homogeneously distributed along the whole circumference for a specific radius. Therefore the maximum axial velocity is lower than in the three dimensional case due to the conservation of mass. There are recirculating currents next to the inlet in the two dimensional cases which increase in magnitude with rotational speed. In the three dimensional case, recirculating flows are located next to the inlet as well, but displaced few degrees in the angle of rotation.

Nevertheless the predicted axial velocity profiles in the two dimensional cases are comparable to the ones obtained in the three dimensional ones in respect to the radial distribution of the flow. The inlet jet remains until the water leaves the centrifuge via the circular weir, so that no plug-shaped profile - resulting in a high active volume -, as observed in the experiments, is formed.

### **4.4 Calculation of the flow pattern (k-**e **RNG model, VOF method, transient solution) overflow setup**

The overflow setup was simulated by using the previously described Volume of Fluid method. The centrifuge is filled with water at the beginning of the simulation. As an initial condition, the water is already spinning with the rigid body rotation. This assumption is justifiable because in the experiments, the centrifuge was slowly filled with water and ran for at least 30 minutes in steady state prior to the measurements. The initial condition reduces the necessary time until the steady state is reached significantly. Figure 12 shows

Fig. 12. Profiles of water-air interface for a throughput of 1 l/min: a) 1000 rpm; b) 9000 rpm

the distribution of the two phases for 1000 (a) und 9000 rpm (b) for a throughput of 1 l/min. Due to the gravitational and centrifugal force, a conical shape of the water surface is developed at 1000 rpm. The centrifugal number outbalances the gravitational force with

Computational Fluid Dynamics (CFD) and

and 5.8 l/min

Discrete Element Method (DEM) Applied to Centrifuges 117

does not widen so that a short-circuit flow towards the outlet is formed, as shown in the contour plot on the left side of Figure 14. There is no significant expansion of the flow profile neither in angular nor in radial direction. The peak of the axial velocity shifts from a radius of 42 mm close to the inlet towards a radius of 36 mm. This behaviour was also predicted in the single-phase model. Nevertheless the simulation misrepresents the true

In order to investigate how the acceleration of the feed affects the flow pattern, a centrifuge with a simpler geometry and a low performance feed accelerator was simulated. The main differences when compared to the tubular bowl centrifuge are the lower length-to-diameter ratio of 0.8 and the feed accelerator. The feed enters a rotating disc and leaves it radially with a certain angular velocity. The water exits the centrifuge through eight boreholes distributed along the circumference at the top of the bowl. Figure 15 shows a scheme of the centrifuge geometry on the left and the results of the volume fraction of air on the right. Water is fed axially through the inlet to the accelerator; there, the water changes its direction and gains in tangential velocity. The transfer of rotational movement occurs just at the plate surfaces, defined as no-slip boundaries, via momentum diffusion. Then, the water exits the accelerator as streams that travel through the air core and enter the rotating liquid pool. The radius of the interface changes 35 % along the height of the centrifuge from the inlet to the

Fig. 15. Left: geometry of the centrifuge; right: contours of volume fraction of air at 2550 rpm

The main flow occurs in the direction of the rotation of the bowl. Both, the liquid pool and the air core rotate. The tangential velocity of the water jet from the inlet accelerator (height 0.0275 m) is lower than the tangential velocity of the bowl. This causes a drag of the whole rotating liquid pool regarding the angular velocity of the bowl. The values of the tangential velocity of the liquid averaged over the circumference are plotted for different heights in Figure 16. The averaged velocity of the liquid in this geometry is approximately 50 % below the rigid body rotation. The values at the wall have to match the rigid body rotation of the bowl to fulfil the no-slip boundary condition. The steep decrease of the tangential velocity

case, in which a significant widening of the inlet jet was observed.

upper part of the centrifuge because of the low rotational speed (454 g).

**4.5 Influence of the inlet geometry on the flow pattern** 

increasing rotational speed, thus the water phase approaches the shape of a cylinder, as shown in Figure 12 b).

The tangential velocity profiles for 1000 and 9000 rpm are shown in Figure 13. The diagram on the right side shows the tangential velocity values averaged over the length of the centrifuge for the water phase only. The computed values are compared to the rigid body rotation. The contours of the tangential velocity at 9000 rpm are shown on the left side of Figure 13. The calculated velocities are in good agreement with the analytical values and with the measurements. At 9000 rpm, the tangential velocity of the water exhibits a small lag when compared to the rigid body rotation. The magnitude of the deviation is within the accuracy of the measurements, so that no statement whether this difference is the true case or not, is possible.

Fig. 13. Tangential velocity profiles, water phase, throughput of 1 l/min; left: contours at 9000 rpm; right: for different rotational speeds and different axial positions

Figure 14 shows the axial velocity profiles for 9000 rpm at a constant throughput of 1 l/min. The axial flow pattern is similar to the one predicted in the single-phase model. The inlet jet

Fig. 14. Axial velocity profiles at 9000 rpm and 1 l/min throughput; left: contours; right: different distances from the inlet

increasing rotational speed, thus the water phase approaches the shape of a cylinder, as

The tangential velocity profiles for 1000 and 9000 rpm are shown in Figure 13. The diagram on the right side shows the tangential velocity values averaged over the length of the centrifuge for the water phase only. The computed values are compared to the rigid body rotation. The contours of the tangential velocity at 9000 rpm are shown on the left side of Figure 13. The calculated velocities are in good agreement with the analytical values and with the measurements. At 9000 rpm, the tangential velocity of the water exhibits a small lag when compared to the rigid body rotation. The magnitude of the deviation is within the accuracy of the measurements, so that no statement whether this difference is the true case

Fig. 13. Tangential velocity profiles, water phase, throughput of 1 l/min; left: contours at

Fig. 14. Axial velocity profiles at 9000 rpm and 1 l/min throughput; left: contours; right:

different distances from the inlet

Figure 14 shows the axial velocity profiles for 9000 rpm at a constant throughput of 1 l/min. The axial flow pattern is similar to the one predicted in the single-phase model. The inlet jet

9000 rpm; right: for different rotational speeds and different axial positions

shown in Figure 12 b).

or not, is possible.

does not widen so that a short-circuit flow towards the outlet is formed, as shown in the contour plot on the left side of Figure 14. There is no significant expansion of the flow profile neither in angular nor in radial direction. The peak of the axial velocity shifts from a radius of 42 mm close to the inlet towards a radius of 36 mm. This behaviour was also predicted in the single-phase model. Nevertheless the simulation misrepresents the true case, in which a significant widening of the inlet jet was observed.

### **4.5 Influence of the inlet geometry on the flow pattern**

In order to investigate how the acceleration of the feed affects the flow pattern, a centrifuge with a simpler geometry and a low performance feed accelerator was simulated. The main differences when compared to the tubular bowl centrifuge are the lower length-to-diameter ratio of 0.8 and the feed accelerator. The feed enters a rotating disc and leaves it radially with a certain angular velocity. The water exits the centrifuge through eight boreholes distributed along the circumference at the top of the bowl. Figure 15 shows a scheme of the centrifuge geometry on the left and the results of the volume fraction of air on the right. Water is fed axially through the inlet to the accelerator; there, the water changes its direction and gains in tangential velocity. The transfer of rotational movement occurs just at the plate surfaces, defined as no-slip boundaries, via momentum diffusion. Then, the water exits the accelerator as streams that travel through the air core and enter the rotating liquid pool. The radius of the interface changes 35 % along the height of the centrifuge from the inlet to the upper part of the centrifuge because of the low rotational speed (454 g).

Fig. 15. Left: geometry of the centrifuge; right: contours of volume fraction of air at 2550 rpm and 5.8 l/min

The main flow occurs in the direction of the rotation of the bowl. Both, the liquid pool and the air core rotate. The tangential velocity of the water jet from the inlet accelerator (height 0.0275 m) is lower than the tangential velocity of the bowl. This causes a drag of the whole rotating liquid pool regarding the angular velocity of the bowl. The values of the tangential velocity of the liquid averaged over the circumference are plotted for different heights in Figure 16. The averaged velocity of the liquid in this geometry is approximately 50 % below the rigid body rotation. The values at the wall have to match the rigid body rotation of the bowl to fulfil the no-slip boundary condition. The steep decrease of the tangential velocity

Computational Fluid Dynamics (CFD) and

**5. Mathematical formulation of the particles** 

the Euler-Euler and the Euler-Lagrange methods.

Discrete Element Method (DEM) Applied to Centrifuges 119

There are several multiphase models available for simulation of a particulate phase. A summary of the different simulation methods for multiphase flows in CFD can be found in (van Wachem & Almstedt, 2003). The multiphase models can be divided into two groups;

The Euler-Euler method considers all phases as continuous ones. For each phase, a momentum conservation equation is solved. In the momentum equation of the dispersed phases, the kinetic theory for granular flows (Lun & Savage, 1987) is used. The intensity of the particle velocity fluctuations determines the stresses, viscosity, and pressure of the solid phase. The interaction between phases is considered by means of momentum exchange coefficients. This numerical method has the advantage that it is suitable for high volume concentrations of the disperse phase. The disadvantage is that the properties of each dispersed particle cannot be simulated and a particle size distribution cannot be considered. Moreover, the stability of the calculation depends on the choice of model parameters and often convergence problems appear. For the simulation of sedimentation and thickening processes there are also some Euler-Euler models based on the Kynch's kinematics sedimentation theory (Kynch, 1952), which were reviewed by Bürger in (Bürger & Wendland, 2001). These models are all based on the measurement of fundamental material properties of the suspension to be separated as explained by (Buscall, 1990; Buscall et al., 1987; Landman & White, 1994). These multiphase models are based on mass balances for the fluid and the disperse phase. The fluid flow is not calculated but they consider the relative velocity between fluid and particles. On this basis, these models are unsuitable for studying the particle behaviour in solid bowl centrifuges where the liquid flow plays a crucial role in particle settling. These models usually describe the sedimentation and consolidation in batch processes in simple-geometry centrifuges. An extension of these models is necessary

to describe the sediment build-up in two and three dimensions correctly.

decide which multiphase model and coupling methodology should be used.

The Euler-Lagrange method solves the continuum conservation equations just for the characterization of the continuous phase. It describes the dispersed phase as mass points and determines the particle trajectories by integrating a force balance for each individual particle. This way, instead of a partial differential equation, only an ordinary differential equation must be solved for the dispersed phase. The weakness of the Euler-Lagrange method is that it is only valid for diluted suspensions if a one-way or a two-way coupling method between continuous and disperse phase is applied. With the one-way coupling, the motion of the discrete phase is calculated based on the velocity field of the continuous phase but the particles have no effect on the pressure and velocity field of the continuous phase. By using a two-way coupling, the equations of continuous and discrete phase are calculated alternatively. The momentum conservation equation of the continuous phase has a source term to consider the momentum exchange with the particulate phase. The maximum value for the volume concentration of the disperse phase is between 5% (Schütz et al., 2007) and 10-12% (ANSYS, 2009), above this value the interaction between particles should be considered. Thus there is a third method, four-way coupling, which involves the interaction within the disperse phase: impacts between particles, particle-wall interactions, van der Waals forces, electrostatic forces, etc. Thus the restriction for the volume fraction is not valid anymore. However, if the number of particles increases, the computational effort rises significantly. Usually the characteristics of the problem and the information that should be obtained from the simulation

near the wall was found to be smoother for a 6 times finer grid. However, the computational effort for the multiphase simulation including the particles with the finer grid was unaffordable. The deviation with respect to the rigid body rotation decreases for lower flow rates because a smaller amount of liquid has to be accelerated. The simulation of this geometry shows how important is to accelerate the feed properly to the angular velocity of the liquid pool before entering it in order to reach the angular velocity of the bowl in the rotating liquid and thus the highest centrifugal acceleration possible.

In contrast to the other centrifuge, a layer with high axial velocity at the gas-liquid interface, as expected by the boundary layer model, has been observed in the simulations as depicted in the diagram on the right in Figure 16. Near the impact height of the inlet water jet two whirls with opposed directions are formed. This can be seen at the negative axial velocity values at the interface for the height 0.02 m, just below the inlet, and the positive axial velocities at he interface for the height 0.0275 m, just above the inlet. Upsides the height of the inlet, a homogeneous axial boundary layer is formed at the interface with a recirculation layer at the bowl wall. This layer exhibits a similar averaged thickness δ as the one proposed by the theoretical model of Gosele (Gösele, 1968) but much thicker as the one calculated following the theory of Reuter (Reuter, 1967) as seen in Table 2. Indeed, it occupies the whole liquid pool until a recirculation layer is formed close to the bowl wall. This is due to the lower tangential velocities which cause a small pressure gradient in the radial direction reducing the stratification of the flow and allowing the fluid to move in the axial direction. The velocity values oscillate with the angle, especially near the inlet at 0.02 m and 0.0275 m and the outlet at 0.0975 m, respectively. Thus the standard deviation is higher at these positions.

Fig. 16. Left: tangential velocity versus radial position; right: axial velocity versus radial position at 2550 rpm and 5.8 l/min


Table 2. Thickness of the axial boundary layer δ obtained in the simulations and calculated analytically at 2550 rpm and a throughput of 5.8 l/min

near the wall was found to be smoother for a 6 times finer grid. However, the computational effort for the multiphase simulation including the particles with the finer grid was unaffordable. The deviation with respect to the rigid body rotation decreases for lower flow rates because a smaller amount of liquid has to be accelerated. The simulation of this geometry shows how important is to accelerate the feed properly to the angular velocity of the liquid pool before entering it in order to reach the angular velocity of the bowl in the

In contrast to the other centrifuge, a layer with high axial velocity at the gas-liquid interface, as expected by the boundary layer model, has been observed in the simulations as depicted in the diagram on the right in Figure 16. Near the impact height of the inlet water jet two whirls with opposed directions are formed. This can be seen at the negative axial velocity values at the interface for the height 0.02 m, just below the inlet, and the positive axial velocities at he interface for the height 0.0275 m, just above the inlet. Upsides the height of the inlet, a homogeneous axial boundary layer is formed at the interface with a recirculation layer at the bowl wall. This layer exhibits a similar averaged thickness δ as the one proposed by the theoretical model of Gosele (Gösele, 1968) but much thicker as the one calculated following the theory of Reuter (Reuter, 1967) as seen in Table 2. Indeed, it occupies the whole liquid pool until a recirculation layer is formed close to the bowl wall. This is due to the lower tangential velocities which cause a small pressure gradient in the radial direction reducing the stratification of the flow and allowing the fluid to move in the axial direction. The velocity values oscillate with the angle, especially near the inlet at 0.02 m and 0.0275 m and the outlet at 0.0975 m, respectively. Thus the standard

Fig. 16. Left: tangential velocity versus radial position; right: axial velocity versus radial

**δsimulated (mm) δReuter (mm) δGösele (mm)**  14.1 1.71 15.8 Table 2. Thickness of the axial boundary layer δ obtained in the simulations and calculated

rotating liquid and thus the highest centrifugal acceleration possible.

deviation is higher at these positions.

position at 2550 rpm and 5.8 l/min

analytically at 2550 rpm and a throughput of 5.8 l/min

### **5. Mathematical formulation of the particles**

There are several multiphase models available for simulation of a particulate phase. A summary of the different simulation methods for multiphase flows in CFD can be found in (van Wachem & Almstedt, 2003). The multiphase models can be divided into two groups; the Euler-Euler and the Euler-Lagrange methods.

The Euler-Euler method considers all phases as continuous ones. For each phase, a momentum conservation equation is solved. In the momentum equation of the dispersed phases, the kinetic theory for granular flows (Lun & Savage, 1987) is used. The intensity of the particle velocity fluctuations determines the stresses, viscosity, and pressure of the solid phase. The interaction between phases is considered by means of momentum exchange coefficients. This numerical method has the advantage that it is suitable for high volume concentrations of the disperse phase. The disadvantage is that the properties of each dispersed particle cannot be simulated and a particle size distribution cannot be considered. Moreover, the stability of the calculation depends on the choice of model parameters and often convergence problems appear. For the simulation of sedimentation and thickening processes there are also some Euler-Euler models based on the Kynch's kinematics sedimentation theory (Kynch, 1952), which were reviewed by Bürger in (Bürger & Wendland, 2001). These models are all based on the measurement of fundamental material properties of the suspension to be separated as explained by (Buscall, 1990; Buscall et al., 1987; Landman & White, 1994). These multiphase models are based on mass balances for the fluid and the disperse phase. The fluid flow is not calculated but they consider the relative velocity between fluid and particles. On this basis, these models are unsuitable for studying the particle behaviour in solid bowl centrifuges where the liquid flow plays a crucial role in particle settling. These models usually describe the sedimentation and consolidation in batch processes in simple-geometry centrifuges. An extension of these models is necessary to describe the sediment build-up in two and three dimensions correctly.

The Euler-Lagrange method solves the continuum conservation equations just for the characterization of the continuous phase. It describes the dispersed phase as mass points and determines the particle trajectories by integrating a force balance for each individual particle. This way, instead of a partial differential equation, only an ordinary differential equation must be solved for the dispersed phase. The weakness of the Euler-Lagrange method is that it is only valid for diluted suspensions if a one-way or a two-way coupling method between continuous and disperse phase is applied. With the one-way coupling, the motion of the discrete phase is calculated based on the velocity field of the continuous phase but the particles have no effect on the pressure and velocity field of the continuous phase. By using a two-way coupling, the equations of continuous and discrete phase are calculated alternatively. The momentum conservation equation of the continuous phase has a source term to consider the momentum exchange with the particulate phase. The maximum value for the volume concentration of the disperse phase is between 5% (Schütz et al., 2007) and 10-12% (ANSYS, 2009), above this value the interaction between particles should be considered. Thus there is a third method, four-way coupling, which involves the interaction within the disperse phase: impacts between particles, particle-wall interactions, van der Waals forces, electrostatic forces, etc. Thus the restriction for the volume fraction is not valid anymore. However, if the number of particles increases, the computational effort rises significantly. Usually the characteristics of the problem and the information that should be obtained from the simulation decide which multiphase model and coupling methodology should be used.

Computational Fluid Dynamics (CFD) and

a wall and instead remains in contact with it.

Fig. 18. Coupling between CFD and DEM

limited due to the computational demand.

**5.1 Particles simulation parameters** 

particles by the fluid.

Discrete Element Method (DEM) Applied to Centrifuges 121

(Gondret et al., 2002), Gondret presents the coefficient of restitution as a function of the Stokes number St (Equation 16), which describes the ratio of the dynamic response time of a particle to flow changes and the characteristic flow time. He determined empirically a critical Stokes number of about 10. Below this value, a particle does not bounce when it hits

> <sup>v</sup> <sup>x</sup> 9

In a centrifugal field, the hydrodynamic forces (drag, lift, torque) play a crucial role and must be included in the calculation of the particle trajectories. Therefore the fluid flow and particles movement have to be computed alternatively. A scheme of the coupling between the CFD software FLUENT and the DEM software EDEM can be seen in Figure 18. At first, the flow will be determined with CFD, then for each position in the computational domain, the velocity, pressure, density and viscosity of the flow are transferred to DEM in order to calculate the hydrodynamic forces and include them in the balance. Afterwards, the fluid flow will be computed again considering the influence of the particles as an additional sink term in the momentum equations. This sink term include the force transmitted to the

Some authors have already simulated such complex multiphase flows in other separation equipment such as filters and hydrocyclones (Chu et al., 2009; Chu & Yu, 2008; Dong et al., 2003; Ni et al., 2006). Nevertheless due to the high velocity gradients and complex interactions between the phases, the simulation of multiphase flow and in particular the sediment behaviour in centrifugal field is still an unsolved challenge. The advantage of this approach is that an accurate description of the flow, the particle paths and the sediment behaviour can be obtained. However, the amount of simulated particles and its size is

The calculation in CFD must be run as an unsteady case when coupled with DEM. The DEM time steps are typically smaller than the CFD time steps in order to correctly capture the contact behaviour. The choice of time steps in DEM simulation is of great importance in order to achieve stability in the calculation but does not perform very long calculations, although still representing the real particle behaviour. Therefore it is necessary to analyse

<sup>1</sup> St s (16)

The Discrete Element Method (DEM) (Cundall & Strack, 1979) is a Lagrangian method originally developed for the simulation of bulk solids and then extended to the entire field of particle technology. This method solves the Newton's equations of motion for translation and rotation of the particles. In these equations, impact forces of particle-particle and particle-wall interactions are considered. Other forces such as electrostatic or Brownian forces can also be implemented if necessary. In order to describe the collisions of particles, the soft-sphere approach was chosen. When two particles collide, they actually deform. This deformation is described by an overlap displacement of the particles. A schematic representation of the forces on the particles in the soft-sphere model is depicted in Figure 17. This approach uses a spring-dashpot-slider system to calculate the behaviour of the particles during the contact time. The soft-sphere model allows contacts of a particle with more neighbouring particles during a time step. The net contact force acting on a particle is added in the case of multiple overlapping. This model for contact forces was first developed by Cundall (Cundall & Strack, 1979), who proposed a parallel linear spring-dashpot model for the normal direction and a parallel linear spring-dashpot in a series with a slider for the tangential direction. A comparison of the performance of several soft-sphere models can be found in (Di Renzo & Di Maio, 2005; Stevens & Hrenya, 2005). There, also the Hertz-Mindlin (Mindlin, 1949) non-linear soft-sphere model used in the present work is analysed. This model proposes a normal force as a function of the normal overlap δn, the Young's modulus of the particle material and the particle radii. The damping force in normal direction depends on the restitution coefficient, the normal stiffness, and the normal component of the relative velocity between the particles. The tangential force depends on the tangential overlap δt, on the shear modulus of the particle material, and on the particle radii; while in the damping force in tangential direction the tangential stiffness and the tangential component of the relative velocity between the particles are considered.

Fig. 17. Contact forces between two particles following the soft-sphere model

While the DEM model was originally used in environments without fluid, increasing efforts have been made to extend this model to the mixture of fluid and particles. For that, a coupling of the DEM calculations with the fluid flow simulation is necessary to be able to include the hydrodynamic forces in the particle simulation. An important aspect in the calculation of particles in a fluid is the energy dissipation into the surrounding fluid. In viscous fluids, there are, in contrast to air, two additional effects that lead to energy losses. The drag force leads to lower particle velocities and when a particle is about to hit a surface or another particle, a deceleration occurs. This is due to the drainage of the fluid, which is located in the gap between the two particles or the particle and the surface. In his work

The Discrete Element Method (DEM) (Cundall & Strack, 1979) is a Lagrangian method originally developed for the simulation of bulk solids and then extended to the entire field of particle technology. This method solves the Newton's equations of motion for translation and rotation of the particles. In these equations, impact forces of particle-particle and particle-wall interactions are considered. Other forces such as electrostatic or Brownian forces can also be implemented if necessary. In order to describe the collisions of particles, the soft-sphere approach was chosen. When two particles collide, they actually deform. This deformation is described by an overlap displacement of the particles. A schematic representation of the forces on the particles in the soft-sphere model is depicted in Figure 17. This approach uses a spring-dashpot-slider system to calculate the behaviour of the particles during the contact time. The soft-sphere model allows contacts of a particle with more neighbouring particles during a time step. The net contact force acting on a particle is added in the case of multiple overlapping. This model for contact forces was first developed by Cundall (Cundall & Strack, 1979), who proposed a parallel linear spring-dashpot model for the normal direction and a parallel linear spring-dashpot in a series with a slider for the tangential direction. A comparison of the performance of several soft-sphere models can be found in (Di Renzo & Di Maio, 2005; Stevens & Hrenya, 2005). There, also the Hertz-Mindlin (Mindlin, 1949) non-linear soft-sphere model used in the present work is analysed. This model proposes a normal force as a function of the normal overlap δn, the Young's modulus of the particle material and the particle radii. The damping force in normal direction depends on the restitution coefficient, the normal stiffness, and the normal component of the relative velocity between the particles. The tangential force depends on the tangential overlap δt, on the shear modulus of the particle material, and on the particle radii; while in the damping force in tangential direction the tangential stiffness and the tangential

component of the relative velocity between the particles are considered.

Fig. 17. Contact forces between two particles following the soft-sphere model

While the DEM model was originally used in environments without fluid, increasing efforts have been made to extend this model to the mixture of fluid and particles. For that, a coupling of the DEM calculations with the fluid flow simulation is necessary to be able to include the hydrodynamic forces in the particle simulation. An important aspect in the calculation of particles in a fluid is the energy dissipation into the surrounding fluid. In viscous fluids, there are, in contrast to air, two additional effects that lead to energy losses. The drag force leads to lower particle velocities and when a particle is about to hit a surface or another particle, a deceleration occurs. This is due to the drainage of the fluid, which is located in the gap between the two particles or the particle and the surface. In his work (Gondret et al., 2002), Gondret presents the coefficient of restitution as a function of the Stokes number St (Equation 16), which describes the ratio of the dynamic response time of a particle to flow changes and the characteristic flow time. He determined empirically a critical Stokes number of about 10. Below this value, a particle does not bounce when it hits a wall and instead remains in contact with it.

$$\text{St} = \frac{1}{9} \frac{\rho\_s \mathbf{v} \times \mathbf{x}}{\mu} \tag{16}$$

In a centrifugal field, the hydrodynamic forces (drag, lift, torque) play a crucial role and must be included in the calculation of the particle trajectories. Therefore the fluid flow and particles movement have to be computed alternatively. A scheme of the coupling between the CFD software FLUENT and the DEM software EDEM can be seen in Figure 18. At first, the flow will be determined with CFD, then for each position in the computational domain, the velocity, pressure, density and viscosity of the flow are transferred to DEM in order to calculate the hydrodynamic forces and include them in the balance. Afterwards, the fluid flow will be computed again considering the influence of the particles as an additional sink term in the momentum equations. This sink term include the force transmitted to the particles by the fluid.

Fig. 18. Coupling between CFD and DEM

Some authors have already simulated such complex multiphase flows in other separation equipment such as filters and hydrocyclones (Chu et al., 2009; Chu & Yu, 2008; Dong et al., 2003; Ni et al., 2006). Nevertheless due to the high velocity gradients and complex interactions between the phases, the simulation of multiphase flow and in particular the sediment behaviour in centrifugal field is still an unsolved challenge. The advantage of this approach is that an accurate description of the flow, the particle paths and the sediment behaviour can be obtained. However, the amount of simulated particles and its size is limited due to the computational demand.

### **5.1 Particles simulation parameters**

The calculation in CFD must be run as an unsteady case when coupled with DEM. The DEM time steps are typically smaller than the CFD time steps in order to correctly capture the contact behaviour. The choice of time steps in DEM simulation is of great importance in order to achieve stability in the calculation but does not perform very long calculations, although still representing the real particle behaviour. Therefore it is necessary to analyse

Computational Fluid Dynamics (CFD) and

results of this simulation are depicted in Figure 19.

Discrete Element Method (DEM) Applied to Centrifuges 123

sediment build-up it is not possible to use periodic boundaries as was done for the simulation of the fluid. Moreover, we chose this geometry with the VOF model for air and liquid in order to investigate the interaction between the multiphase model and the use of DEM to calculate the particles. The calculations were made with the software ANSYS

The coupled simulation was performed with a rotational speed of 2550 rpm. For this rotational speed, the Stokes settling velocity of the particles is by two orders of magnitude higher than the Stokes settling velocity of the quartz particles used in the experiments. The

Fig. 19. Particles and path lines simulated with EDEM coupled with a two-phase simulation

The particles are injected at the inlet pipe with the same axial velocity as the fluid. They hit the top of the feed accelerator and gain in radial velocity. Afterwards, they move through the air core until they hit the liquid pool. There, under the effect of the centrifugal force, the particles perform spiral path lines, as shown in Figure 19 on the right, until reaching the wall or already settled particles. The lower plot in Figure 19 shows a vertical plane of the centrifuge with the particles in the region near the feed accelerator where the sediment is built. Due to the drag force of the axial flow in the rotating pool, the particles move first to higher axial positions. Then, near the wall, they move backwards because of the recirculation layer. The evolution of the sediment can be observed in Figure 20 as a function

of a centrifuge in FLUENT and coloured by velocity magnitude

FLUENT 12.0.16 and EDEM 2.3, which include a coupling feature for both programs.

the forces acting on the particles to determine the different time scales. In the case presented, the Hertzian contact time - the total time of contact in Hertz theory of elastic collision -, the particle sedimentation time - the time required for a particle to settle one diameter -, and the particle relaxation time - the time for a particle to adapt to the surrounding fluid flow - are considered. On the basis of these, the smallest time step must be chosen so that all physical phenomena are represented. In the simulations presented, the Hertzian contact time is the smallest time scale and, thus a time step tDEM = 10-6 s was chosen in the DEM software, while in the CFD a time step tCFD = 10-4 s was sufficient.

The size of the particles was defined as a logarithmic normal distribution with a mean diameter of 200 µm and a standard deviation of σ = 10 µm. Simulations with smaller particles lead to divergence of the calculation despite of the small time steps set. Thus simulations with lower densities of the particles and rotational speeds were performed to obtain similar settling velocities as with small quartz particles (mean diameter of 2.07 µm, s=2.15 µm) which were used in the experiments. In order to reproduce the experimental concentration of approximately 0.5 vol. %, 100,000 particles per second were injected at the inlet of the centrifuge for the simulation. In future work, the comparison of predicted and experimentally determined sediment build-up will show whether the simulation methodology leads to correct predictions.

For the contact models, a parameter setting is needed to represent the real particle behaviour. The information needed is mainly about the physical properties of the particle system because the spring and damper constants are calculated as a function of the shear modulus and the Poisson's ratio of the particle and wall materials. The coefficient of restitution depends on the Stokes particle number. The friction coefficient is an empirically determinable property, which can be found in the literature for various material combinations. Some authors have already determined this friction coefficient for sand and glass particles for DEM (Asaf et al., 2007; Li et al., 2005). Here, shear experiments of a sediment made of quartz particles were performed in a Jenike shear cell under normal pressure in order to calculate the friction coefficient as was done by (Hartl & Ooi, 2008).


Table 3. Parameters used in the DEM simulation

### **6. Prediction of the sediment build-up by means of DEM**

The geometry of the centrifuge in which the particle trajectories and the sediment build-up were simulated is shown in Figure 15. For the simulation purpose of predicting the

the forces acting on the particles to determine the different time scales. In the case presented, the Hertzian contact time - the total time of contact in Hertz theory of elastic collision -, the particle sedimentation time - the time required for a particle to settle one diameter -, and the particle relaxation time - the time for a particle to adapt to the surrounding fluid flow - are considered. On the basis of these, the smallest time step must be chosen so that all physical phenomena are represented. In the simulations presented, the Hertzian contact time is the smallest time scale and, thus a time step tDEM = 10-6 s was chosen in the DEM software,

The size of the particles was defined as a logarithmic normal distribution with a mean diameter of 200 µm and a standard deviation of σ = 10 µm. Simulations with smaller particles lead to divergence of the calculation despite of the small time steps set. Thus simulations with lower densities of the particles and rotational speeds were performed to obtain similar settling velocities as with small quartz particles (mean diameter of 2.07 µm, s=2.15 µm) which were used in the experiments. In order to reproduce the experimental concentration of approximately 0.5 vol. %, 100,000 particles per second were injected at the inlet of the centrifuge for the simulation. In future work, the comparison of predicted and experimentally determined sediment build-up will show whether the simulation

For the contact models, a parameter setting is needed to represent the real particle behaviour. The information needed is mainly about the physical properties of the particle system because the spring and damper constants are calculated as a function of the shear modulus and the Poisson's ratio of the particle and wall materials. The coefficient of restitution depends on the Stokes particle number. The friction coefficient is an empirically determinable property, which can be found in the literature for various material combinations. Some authors have already determined this friction coefficient for sand and glass particles for DEM (Asaf et al., 2007; Li et al., 2005). Here, shear experiments of a sediment made of quartz particles were performed in a Jenike shear cell under normal pressure in order to calculate the friction coefficient as was done by (Hartl & Ooi, 2008).

> **Parameter Value**  Particle radius (µm) 200 Particle density (kg/m3) 1016 Shear modulus particle (Pa) 2.2·108 Shear modulus wall (Pa) 7.0·1010 Poisson's ratio particle (-) 0.25 Poisson's ratio wall (-) 0.30 Restitution coefficient (-) 0.010 Static friction coefficient (-) 0.787 Rolling friction coefficient (-) 0.100

The geometry of the centrifuge in which the particle trajectories and the sediment build-up were simulated is shown in Figure 15. For the simulation purpose of predicting the

while in the CFD a time step tCFD = 10-4 s was sufficient.

methodology leads to correct predictions.

Table 3. Parameters used in the DEM simulation

**6. Prediction of the sediment build-up by means of DEM** 

sediment build-up it is not possible to use periodic boundaries as was done for the simulation of the fluid. Moreover, we chose this geometry with the VOF model for air and liquid in order to investigate the interaction between the multiphase model and the use of DEM to calculate the particles. The calculations were made with the software ANSYS FLUENT 12.0.16 and EDEM 2.3, which include a coupling feature for both programs.

The coupled simulation was performed with a rotational speed of 2550 rpm. For this rotational speed, the Stokes settling velocity of the particles is by two orders of magnitude higher than the Stokes settling velocity of the quartz particles used in the experiments. The results of this simulation are depicted in Figure 19.

Fig. 19. Particles and path lines simulated with EDEM coupled with a two-phase simulation of a centrifuge in FLUENT and coloured by velocity magnitude

The particles are injected at the inlet pipe with the same axial velocity as the fluid. They hit the top of the feed accelerator and gain in radial velocity. Afterwards, they move through the air core until they hit the liquid pool. There, under the effect of the centrifugal force, the particles perform spiral path lines, as shown in Figure 19 on the right, until reaching the wall or already settled particles. The lower plot in Figure 19 shows a vertical plane of the centrifuge with the particles in the region near the feed accelerator where the sediment is built. Due to the drag force of the axial flow in the rotating pool, the particles move first to higher axial positions. Then, near the wall, they move backwards because of the recirculation layer. The evolution of the sediment can be observed in Figure 20 as a function

Computational Fluid Dynamics (CFD) and

taken from the simulation domain.

Discrete Element Method (DEM) Applied to Centrifuges 125

position calculation of the particles. Thus in Figure 21 on the left, the particles calculated in FLUENT acquire higher axial positions when compared to the particles calculated with EDEM. The reason for this disagreement is the difference in the axial flow pattern, where the velocities oscillate until a quasi-steady state is achieved. However, qualitatively both particle trajectories are in good agreement. The tangential velocity of the particles is represented in Figure 21 on the right. The tangential velocity of the particles is almost zero in the regions where there is no liquid. Then, they are accelerated in tangential direction until reaching the bowl wall. There, a deviation between FLUENT and EDEM appears. The particles in EDEM hit the sediment and the tangential velocity diminishes because the wall rotation is not taken into account here, but in FLUENT, no particle accumulation is considered and, thus all particles are calculated until they reach the wall. Once they have reached the wall, they are

The location and quantity of the particles are stored as a variable in FLUENT representing the volume fraction of the particles, as can be seen in Figure 22 on the left. However, for the flow calculations they are considered as points where a momentum exchange occurs. The drag force of the flow acting on the particles is subtracted from the momentum equation diminishing the flow in regions where the particles accumulate. Figure 22 on the right represents the volume fraction of air in the centrifuge. The effect of the particles on the flow can be seen as the water jet coming from the inlet accelerator falls down in the axial position until it reaches the liquid pool due to the particles weight and the momentum sink term.

Fig. 21. Comparison of the particle trajectories obtained in FLUENT and EDEM: left: axial position versus radius; right: tangential velocity of the particles and the fluid versus radius

of the calculated particles for two different rotational speeds of the liquid pool. For the case with lower rotational speed, 686 rpm, the first settled particles deposit at the wall and are pushed to lower axial positions because of the recirculation layer near the wall and due to the lower normal pressure, 3 kPa, acting on the sediment, which allows the influence of the gravity. Then, particles deposit on the sediment and it growths from the wall to the inner part of the centrifuge in radial direction. The sliding of the sediment in this case is due to the relationship of low normal pressure and high shear stress acting on it. In the case of higher rotational speed in the liquid pool, 1937 rpm, the particles travel radially direct to the bowl wall and the sediment formed is flat. The particles settling in the next time steps slide to the side of the sediment which grows along the wall rather than in the radial direction.

In this study the value of the friction coefficient between particles and between particles and wall was constant. Changes of this parameter would also lead to changes in the sediment behaviour under certain conditions of normal pressure and shear stress. A further investigation in this field, including the experimental determination of the friction coefficient for different materials and particle size distributions, is needed. Here only the basis of the coupling methodology was studied.

Fig. 20. Sediment built at the bowl wall as a function of the particles simulated with EDEM coupled with a two-phase simulation of a centrifuge in FLUENT and coloured by velocity magnitude for two different rotational speeds

In order to ensure that the coupling between both programs works and that the hydrodynamic forces are correctly passed from FLUENT to EDEM, the particle trajectories calculated with the coupling and with FLUENT were evaluated. Although the time steps and the mathematical method used in FLUENT to calculate the particles is different from the time steps and algorithm used in EDEM, the results are comparable as shown in Figure 21. Another different is that in FLUENT, the particles were calculated using a static flow field at a certain time step whereas in EDEM, the flow field was calculated alternating to the

of the calculated particles for two different rotational speeds of the liquid pool. For the case with lower rotational speed, 686 rpm, the first settled particles deposit at the wall and are pushed to lower axial positions because of the recirculation layer near the wall and due to the lower normal pressure, 3 kPa, acting on the sediment, which allows the influence of the gravity. Then, particles deposit on the sediment and it growths from the wall to the inner part of the centrifuge in radial direction. The sliding of the sediment in this case is due to the relationship of low normal pressure and high shear stress acting on it. In the case of higher rotational speed in the liquid pool, 1937 rpm, the particles travel radially direct to the bowl wall and the sediment formed is flat. The particles settling in the next time steps slide to the

side of the sediment which grows along the wall rather than in the radial direction.

basis of the coupling methodology was studied.

magnitude for two different rotational speeds

In this study the value of the friction coefficient between particles and between particles and wall was constant. Changes of this parameter would also lead to changes in the sediment behaviour under certain conditions of normal pressure and shear stress. A further investigation in this field, including the experimental determination of the friction coefficient for different materials and particle size distributions, is needed. Here only the

Fig. 20. Sediment built at the bowl wall as a function of the particles simulated with EDEM coupled with a two-phase simulation of a centrifuge in FLUENT and coloured by velocity

In order to ensure that the coupling between both programs works and that the hydrodynamic forces are correctly passed from FLUENT to EDEM, the particle trajectories calculated with the coupling and with FLUENT were evaluated. Although the time steps and the mathematical method used in FLUENT to calculate the particles is different from the time steps and algorithm used in EDEM, the results are comparable as shown in Figure 21. Another different is that in FLUENT, the particles were calculated using a static flow field at a certain time step whereas in EDEM, the flow field was calculated alternating to the position calculation of the particles. Thus in Figure 21 on the left, the particles calculated in FLUENT acquire higher axial positions when compared to the particles calculated with EDEM. The reason for this disagreement is the difference in the axial flow pattern, where the velocities oscillate until a quasi-steady state is achieved. However, qualitatively both particle trajectories are in good agreement. The tangential velocity of the particles is represented in Figure 21 on the right. The tangential velocity of the particles is almost zero in the regions where there is no liquid. Then, they are accelerated in tangential direction until reaching the bowl wall. There, a deviation between FLUENT and EDEM appears. The particles in EDEM hit the sediment and the tangential velocity diminishes because the wall rotation is not taken into account here, but in FLUENT, no particle accumulation is considered and, thus all particles are calculated until they reach the wall. Once they have reached the wall, they are taken from the simulation domain.

The location and quantity of the particles are stored as a variable in FLUENT representing the volume fraction of the particles, as can be seen in Figure 22 on the left. However, for the flow calculations they are considered as points where a momentum exchange occurs. The drag force of the flow acting on the particles is subtracted from the momentum equation diminishing the flow in regions where the particles accumulate. Figure 22 on the right represents the volume fraction of air in the centrifuge. The effect of the particles on the flow can be seen as the water jet coming from the inlet accelerator falls down in the axial position until it reaches the liquid pool due to the particles weight and the momentum sink term.

Fig. 21. Comparison of the particle trajectories obtained in FLUENT and EDEM: left: axial position versus radius; right: tangential velocity of the particles and the fluid versus radius

Computational Fluid Dynamics (CFD) and

Discrete Element Method (DEM) Applied to Centrifuges 127

**7. Evaluation of the computational demand for predicting flow in centrifuges**  The prediction of flow pattern in centrifuges requires a high computational demand that increases with the maximum velocity in the system. Accordingly to Equation 17, the highest

velocity is reached at the bowl wall and increases linearly with the rotational speed.

2 60 *n v r* 

Figure 24 shows the computational demand in hsim/(hflow time 3.6 10-3) (dimensionless) against circumferential velocity in m/s for different simulations which include two and three dimensional problems, one and two-phase simulations, and different models to predict turbulent flow. The simulations were performed with Win xp 32 bit desktop computers using an Intel Core 2 Duo E4500 2.2 GHz processor, 2 GB DDR 2 memory running at 266 MHz (2 x 1GB Kingston) and a ASrock motherboard (dual channel memory). All data are averaged values. For reasons of readability, the standard deviation is not stated, because it is too small to be included in the plot. Evaluating the diagram on the left side of Figure 24, it is remarkable how costly the VOF method for simulating two-phase flow for the two and three dimensional cases is. Moving from one to two phases causes the same increase in computational demand as moving from a two to a three dimensional grid and simultaneously increasing the cell number by one order of magnitude. However, the size of the grid does not necessarily increase the computational time. The convergence in a fine grid is faster and the influence of the rotating geometry overcomes the one of the grid when

Fig. 24. Correlation of computational demand in hsim/(hflow time 3.6 10-3) and circumferential velocity: for one and two-phase modelling, two and three-dimensional geometries, and different sizes of the geometry (left); dependence on turbulent model used (right), 715000

cells, three-dimensional and single-phase model

(17)

Fig. 22. Contours of the volume fraction in the FLUENT simulation domain; left: particles; right: air

The effect of the sediment on the flow pattern is shown in Figure 23, where the tangential and axial velocity values are represented. The tangential velocity values are lower than in the simulation without particles because of the momentum exchange between the fluid and the particles. Furthermore the tangential velocity is lower in the lower part of the centrifuge, where most of the particles can be found. The interface radius is not constant over the height of the centrifuge, as can be observed in Figure 22 on the right. Therefore the values of the tangential and axial velocity of the water in Figure 23 begin for lower radii in the lower part of the centrifuge. The form of the axial flow pattern does not change qualitatively except for the height where the particles accumulate, around 0.02 m. There, the axial velocity takes approximately the value zero due to the exchange of momentum between fluid and particles. At the wall, where the sediment is located, there is no axial flow in contrast to the simulation without particles. This fact shows how important is to consider the sediment build-up in the multiphase flow simulation of a centrifuge. However this also carries an increase of the computational effort.

Fig. 23. Averaged values over the circumferential position for different axial positions of the centrifuge obtained in FLUENT coupled with EDEM, left: tangential velocity values versus radius; right: axial velocity versus radius

Fig. 22. Contours of the volume fraction in the FLUENT simulation domain; left: particles;

The effect of the sediment on the flow pattern is shown in Figure 23, where the tangential and axial velocity values are represented. The tangential velocity values are lower than in the simulation without particles because of the momentum exchange between the fluid and the particles. Furthermore the tangential velocity is lower in the lower part of the centrifuge, where most of the particles can be found. The interface radius is not constant over the height of the centrifuge, as can be observed in Figure 22 on the right. Therefore the values of the tangential and axial velocity of the water in Figure 23 begin for lower radii in the lower part of the centrifuge. The form of the axial flow pattern does not change qualitatively except for the height where the particles accumulate, around 0.02 m. There, the axial velocity takes approximately the value zero due to the exchange of momentum between fluid and particles. At the wall, where the sediment is located, there is no axial flow in contrast to the simulation without particles. This fact shows how important is to consider the sediment build-up in the multiphase flow simulation of a centrifuge. However this also carries an

Fig. 23. Averaged values over the circumferential position for different axial positions of the centrifuge obtained in FLUENT coupled with EDEM, left: tangential velocity values versus

right: air

increase of the computational effort.

radius; right: axial velocity versus radius

#### **7. Evaluation of the computational demand for predicting flow in centrifuges**

The prediction of flow pattern in centrifuges requires a high computational demand that increases with the maximum velocity in the system. Accordingly to Equation 17, the highest velocity is reached at the bowl wall and increases linearly with the rotational speed.

$$v = \frac{2 \cdot \pi \cdot n}{60} \cdot r \tag{17}$$

Figure 24 shows the computational demand in hsim/(hflow time 3.6 10-3) (dimensionless) against circumferential velocity in m/s for different simulations which include two and three dimensional problems, one and two-phase simulations, and different models to predict turbulent flow. The simulations were performed with Win xp 32 bit desktop computers using an Intel Core 2 Duo E4500 2.2 GHz processor, 2 GB DDR 2 memory running at 266 MHz (2 x 1GB Kingston) and a ASrock motherboard (dual channel memory). All data are averaged values. For reasons of readability, the standard deviation is not stated, because it is too small to be included in the plot. Evaluating the diagram on the left side of Figure 24, it is remarkable how costly the VOF method for simulating two-phase flow for the two and three dimensional cases is. Moving from one to two phases causes the same increase in computational demand as moving from a two to a three dimensional grid and simultaneously increasing the cell number by one order of magnitude. However, the size of the grid does not necessarily increase the computational time. The convergence in a fine grid is faster and the influence of the rotating geometry overcomes the one of the grid when

Fig. 24. Correlation of computational demand in hsim/(hflow time 3.6 10-3) and circumferential velocity: for one and two-phase modelling, two and three-dimensional geometries, and different sizes of the geometry (left); dependence on turbulent model used (right), 715000 cells, three-dimensional and single-phase model

Computational Fluid Dynamics (CFD) and

steps according to the controlling forces and phenomena.

data, especially in terms of the axial velocity profile.

C2 resistance coefficient, m-1

volumetric force, N m-3

tangential force, N

 p pressure, Pa (N m-2) pc permeability, m2 rbowl bowl radius, m rweir weir radius, m

drag force, N

normal force, N

dV difference in mass balance, %

k turbulent kinetic energy, m2 s-2

n revolutions per minute, rpm

rboundary layer radius of boundary layer, m

L length of the rotor of the centrifuge, m

Cμ constant of the k-ε turbulence model, dimensionless

cycle, can be used.

**9. Notation** 

F 

Fd

Ft 

Fn

Discrete Element Method (DEM) Applied to Centrifuges 129

it proves that the coupling method can be applied for rotating geometries such as centrifuges. The particle trajectories and the sediment build-up agree qualitatively with the expected results. Particles travel through the air core and enter the liquid pool, where they describe spiral paths until they reach the wall or already settled particles to form a sediment. The volume occupied by the particles has not been considered in the calculation of the fluid but the effect of the particles is taken into account by a sink term in the momentum equation. Thus in the flow pattern obtained in FLUENT, the fluid velocity is approximately zero where the sediment is located. The consideration of the particle volume in the calculation of the fluid is a gap in the methodology. A simulation with three Euler phases air, liquid and particles - calculated in EDEM could be a possible but computationally very expensive solution. The method is complex and there is still some effort needed to reduce the particle size without significantly increasing computing time. A possibility to reduce the computational time is to divide the computational domain into regions with different time

There is a substantial increase in computational demand with the increase in rotational speed. This limits the predictions to low rotational speeds when standard workstation computers are used for the calculations. If supercomputers are accessible as it will be the case in research at universities, predictions for high rotational speed are possible. But before advancing to higher rotational speeds it is important to improve the accuracy of the simulations of flow in centrifugal force fields and validate the models with experimental

In order to run a sanity check for the simulations in FLUENT or EDEM, it is not sufficient to evaluate the residuals and mass balance. For new methodologies or boundary conditions, for example when simulating centrifuges, it is necessary to compare different turbulent models, transient and stationary calculation methods and most important to validate the simulation methodology proposed with experimental data. Once validated, geometry and parameters can be changed and the advantage of CFD, which is a fast change to evaluation

simulating centrifuges. This is concluded by comparing the three dimensional, single-phase calculations with a grid size of 215000 cells with the simulations with a grid of 715000 cells. The same effect is observed when comparing the three dimensional, two-phase simulations with different cell sizes and different size of geometry. All data is well fitted with a single function for each of the two compared sets of calculations. Looking at the differences between the turbulent models, all simulations exhibit a comparable computational demand. The RSM is the model with the highest computational demand whereas the large-eddy simulation needs less computational resources than all other models evaluated. This is unexpected since the LES is usually rated as more demanding than the RSM. This is explainable due to the single equation modelling of turbulence in the subgrid, the low degree of turbulence in the case presented and the fine mesh that was constant in all cases presented in the diagram on the right side of Figure 24. The k-e RNG and the k-w simulations require comparable computational resources. This is in agreement with the theory, because both models are similar to each other except of the modelling of the near wall region.

### **8. Conclusions**

The results presented show that the prediction of the flow pattern in centrifuges by means of computational fluid dynamics give good results for evaluation of the tangential acceleration efficiency. The predictions are in good agreement with the measurements conducted with Laser Doppler Anemometry. Nevertheless the axial fluid profiles do not depict the true case accurately. The calculated axial velocity profiles give reasonable results close behind the inlet, but the widening of the inlet jet is not distinguishable in the simulations. The axial velocity profile is essential for the calculation of the cut size and the sediment distribution in the centrifuge, thus any significant deviation from the true case is not acceptable. It was demonstrated that the Reynolds Stress Model (RSM) and the large-eddy simulation (LES) do not increase the accuracy of the predictions when compared to the results obtained using the k-e RNG and k-w model, but increase the computational demand in case of the RSM. All models predict a short-circuit flow between inlet and outlet. The flow pattern does not change with increasing rotational speed or throughputs, only the maximum fluid velocity increases with higher throughputs as expected. This behaviour was observed in single and two-phase simulations.

The geometry of the feed accelerator is decisive for the axial flow pattern developed in the centrifuge. A radial feed impacting on the surface of the liquid pool leads to an axial boundary layer flow similar to the one predicted and experimentally verified by other authors such as Glinka (Glinka, 1983). The design of the feed accelerator and its efficiency is essential for a good centrifugal performance. If the liquid is not accelerated to the angular velocity of the bowl, the rotating pool exhibits a significant lag in tangential velocity when compared to the rigid body rotation. The proposed feeding system via rotating inlet bores leads to an efficient tangential acceleration of the fed suspension.

Particles with a diameter of 200 µm were calculated for the first time in a centrifuge using a coupling between the simulation of the flow in FLUENT and the simulation of the particles in EDEM. The results show a good agreement between particles simulated in FLUENT only considering the hydrodynamic forces and the results from EDEM, which implies that the hydrodynamic forces are correctly transferred from one program to the other. Furthermore it proves that the coupling method can be applied for rotating geometries such as centrifuges. The particle trajectories and the sediment build-up agree qualitatively with the expected results. Particles travel through the air core and enter the liquid pool, where they describe spiral paths until they reach the wall or already settled particles to form a sediment. The volume occupied by the particles has not been considered in the calculation of the fluid but the effect of the particles is taken into account by a sink term in the momentum equation. Thus in the flow pattern obtained in FLUENT, the fluid velocity is approximately zero where the sediment is located. The consideration of the particle volume in the calculation of the fluid is a gap in the methodology. A simulation with three Euler phases air, liquid and particles - calculated in EDEM could be a possible but computationally very expensive solution. The method is complex and there is still some effort needed to reduce the particle size without significantly increasing computing time. A possibility to reduce the computational time is to divide the computational domain into regions with different time steps according to the controlling forces and phenomena.

There is a substantial increase in computational demand with the increase in rotational speed. This limits the predictions to low rotational speeds when standard workstation computers are used for the calculations. If supercomputers are accessible as it will be the case in research at universities, predictions for high rotational speed are possible. But before advancing to higher rotational speeds it is important to improve the accuracy of the simulations of flow in centrifugal force fields and validate the models with experimental data, especially in terms of the axial velocity profile.

In order to run a sanity check for the simulations in FLUENT or EDEM, it is not sufficient to evaluate the residuals and mass balance. For new methodologies or boundary conditions, for example when simulating centrifuges, it is necessary to compare different turbulent models, transient and stationary calculation methods and most important to validate the simulation methodology proposed with experimental data. Once validated, geometry and parameters can be changed and the advantage of CFD, which is a fast change to evaluation cycle, can be used.

### **9. Notation**

128 Applied Computational Fluid Dynamics

simulating centrifuges. This is concluded by comparing the three dimensional, single-phase calculations with a grid size of 215000 cells with the simulations with a grid of 715000 cells. The same effect is observed when comparing the three dimensional, two-phase simulations with different cell sizes and different size of geometry. All data is well fitted with a single function for each of the two compared sets of calculations. Looking at the differences between the turbulent models, all simulations exhibit a comparable computational demand. The RSM is the model with the highest computational demand whereas the large-eddy simulation needs less computational resources than all other models evaluated. This is unexpected since the LES is usually rated as more demanding than the RSM. This is explainable due to the single equation modelling of turbulence in the subgrid, the low degree of turbulence in the case presented and the fine mesh that was constant in all cases presented in the diagram on the right side of Figure 24. The k-e RNG and the k-w simulations require comparable computational resources. This is in agreement with the theory, because both models are

The results presented show that the prediction of the flow pattern in centrifuges by means of computational fluid dynamics give good results for evaluation of the tangential acceleration efficiency. The predictions are in good agreement with the measurements conducted with Laser Doppler Anemometry. Nevertheless the axial fluid profiles do not depict the true case accurately. The calculated axial velocity profiles give reasonable results close behind the inlet, but the widening of the inlet jet is not distinguishable in the simulations. The axial velocity profile is essential for the calculation of the cut size and the sediment distribution in the centrifuge, thus any significant deviation from the true case is not acceptable. It was demonstrated that the Reynolds Stress Model (RSM) and the large-eddy simulation (LES) do not increase the accuracy of the predictions when compared to the results obtained using the k-e RNG and k-w model, but increase the computational demand in case of the RSM. All models predict a short-circuit flow between inlet and outlet. The flow pattern does not change with increasing rotational speed or throughputs, only the maximum fluid velocity increases with higher throughputs as expected. This behaviour was observed in single and

The geometry of the feed accelerator is decisive for the axial flow pattern developed in the centrifuge. A radial feed impacting on the surface of the liquid pool leads to an axial boundary layer flow similar to the one predicted and experimentally verified by other authors such as Glinka (Glinka, 1983). The design of the feed accelerator and its efficiency is essential for a good centrifugal performance. If the liquid is not accelerated to the angular velocity of the bowl, the rotating pool exhibits a significant lag in tangential velocity when compared to the rigid body rotation. The proposed feeding system via rotating inlet bores

Particles with a diameter of 200 µm were calculated for the first time in a centrifuge using a coupling between the simulation of the flow in FLUENT and the simulation of the particles in EDEM. The results show a good agreement between particles simulated in FLUENT only considering the hydrodynamic forces and the results from EDEM, which implies that the hydrodynamic forces are correctly transferred from one program to the other. Furthermore

similar to each other except of the modelling of the near wall region.

leads to an efficient tangential acceleration of the fed suspension.

**8. Conclusions** 

two-phase simulations.


Computational Fluid Dynamics (CFD) and

*Ansys Fluent 12.*,

321-340

2-3, 101-116

*Mechanics*, 24, 2, 183-202

*Technology*, 179, 3, 104-114

*Geotechnique*, 29, 1, 47-65

*Mashinostroenie*, 6, 13-14

*Physics of Fluids*, 14, 2, 643-652

1312

657-659

Discrete Element Method (DEM) Applied to Centrifuges 131

ANSYS. (2006). Chapter 23.3.2 Volume Fraction Equation, In: *User's Guide Ansys Fluent 6.3* ANSYS. (2009). Chapter 4.12.2 Standard Wall Functions, In: *Theory Guide Ansys Fluent 12.0* ANSYS. (2009). Chapter 18.4.1 Discretization of the Momentum Equation, In: *Theory Guide* 

ANSYS. (2009). Chapter 23. Modeling Discrete Phase, In: *User's Guide Ansys Fluent 12.0* Asaf, Z., Rubinstein, D. & Shmulevich, I. (2007). Determination of discrete element model parameters required for soil tillage. *Soil & Tillage Research*, 92, 1-2, 227-242 Bass, E. (1959). Strömungen im Fliehkraftfeld I. *Periodica Polytechics chem. Ingenieurwes*, 3,

Bass, E. (1959). Strömungen im Fliehkraftfeld II - Absetzsicherheit von Röhrenzentrifugen.

Brantley, J. N., Willis, D. D., Breillatt, J. P., Gibson, R. F., Patrick, L. C. & Anderson, N. G.

Brennan, M. S., M., N. & Holtham, P. N. (2007). Multiphase modelling of hydrocyclones –

Bürger, R. & Wendland, W. L. (2001). Sedimentation and suspension flows: Historical

Buscall, R. (1990). The sedimentation of concentrated colloidal suspensions. *Colloids and Surfaces Papers Presented at the Society of Chemical Industry Meeting*, 43, 1, 33-53 Buscall, R., Mills, P. D. A., Stewart, R. F., Sutton, D., White, L. R. & Yates, G. E. (1987). The

Chu, K. W., Wang, B., Yu, A. B. & Vince, A. (2009). CFD-DEM modelling of multiphase flow

Chu, K. W. & Yu, A. B. (2008). Numerical simulation of complex particle-fluid flows. *Powder* 

Cundall, P. A. & Strack, O. D. L. (1979). Discrete Numerical-Model for Granular Assemblies.

Di Renzo, A. & Di Maio, F. P. (2005). An improved integral non-linear model for the contact

Glinka, U. (1983). Die Strömung in Überlaufzentrifugen - Neue Ergebnisse mit einem elektrolytischen Markierungsverfahren. *Verfahrenstechnik*, 17, 5, 315-318 Golovko, Y. D. (1969). Liquid Lag in Tubular Centrifuge Rotors. *Khimicheskoe i Neftyanoe* 

Gondret, P., Lance, M. & Petit, L. (2002). Bouncing motion of spherical particles in fluids.

Gösele, W. (1968). Schichtströmung in Röhrenzentrifugen. *Chemie IngenieurTechnik*, 40, 13,

of particles in distinct element simulations. *Chemical Engineering Science*, 60, 5, 1303-

in dense medium cyclones. *Powder Technology*, 193, 3, 235-247

(1970). K-series centrifuges : IV. Measurement and control of temperature.

perspective and some recent developments. *Journal of Engineering Mathematics*, 41,

rheology of strongly-flocculated suspensions. *Journal of Non-Newtonian Fluid* 

for vaccine purification. *Analytical Biochemistry*, 32, 3, 460-494

Bass, E. (1962). Strömungs- und Absetzvorgänge in Röhrenzentrifugen. 249-264

prediction of cut-size. *Minerals Engineering*, 20, 395-406

*Periodica Polytechics chem. Ingenieurwes*, 4, 41-61

*Analytical Biochemistry*, 36, 2, 434-442

Development of the K-II continuous-sample-flow-with-banding centrifuge system


### **10. Acknowledgments**

We acknowledge support by Deutsche Forschungsgemeinschaft and Open Access Publishing Fund of Karlsruhe Institute of Technology.

### **11. References**

Anderson, N. G., Waters, D. A., Nunley, C. E., Gibson, R. F., Schilling, R. M., Denny, E. C., Cline, G. B., Babelay, E. F. & Perardi, T. E. (1969). K-Series Centrifuges I.

Si sink term, N m-3

vax axial velocity, m s-1

δn normal overlap, m δt tangential overlap, m

 density, kg m-3 s solids density, kg m-3

velocity, m s-1

 u velocity in x direction, m s-1 v velocity in y direction, m s-1

 w velocity in z direction, m s-1 x particle diameter, m xcell cell length, m

α volume fraction, dimensionless

εs sediment porosity, dimensionless µ dynamic viscosity, kg m-1 s-1 µt turbulent viscosity, kg m-1 s-1

τt turbulent stress tensor, Pa (N m-2) τw wall shear stress, Pa (N m-2) i Poissons ratio, dimensionless

' fluctuating component of a variable

Publishing Fund of Karlsruhe Institute of Technology.

 LES large-eddy simulation RBR rigid body rotation RSM Reynolds Stress Model VOF Volume of Fluid St Stokes number

δ thickness of the axial boundary layer, mm

ε turbulent kinetic energy dissipation rate, m2 s-3

density difference between liquid and solid, kg m-3

specific turbulent kinetic energy dissipation rate, s-1

We acknowledge support by Deutsche Forschungsgemeinschaft and Open Access

Anderson, N. G., Waters, D. A., Nunley, C. E., Gibson, R. F., Schilling, R. M., Denny, E. C.,

Cline, G. B., Babelay, E. F. & Perardi, T. E. (1969). K-Series Centrifuges I.

y+ characteristic width of viscous sublayer, dimensionless

t time, s

v 

Greek letters

Indices

q phase

**10. Acknowledgments** 

**11. References** 

Abbreviations

Development of the K-II continuous-sample-flow-with-banding centrifuge system for vaccine purification. *Analytical Biochemistry*, 32, 3, 460-494


Computational Fluid Dynamics (CFD) and

Discrete Element Method (DEM) Applied to Centrifuges 133

Reuter, H. (1967). Sedimentation in der Überlaufzentrifuge. *Chemie-Ing. Tech.*, 39, 9, 548-553 Romani Fernandez, X. & Nirschl, H. (2009). Multiphase CFD Simulation of a Solid Bowl

Saunders, E. (1948). Nomograph for Particle Size Determination with the Sharples

Schaeflinger, U. & Stibi, H. (1991). Centrifugal separation of a mixture in a rotating bucket II.

Schaeflinger, U., Stibi, H. (1991). Centrifugal separation of a mixture in a rotating bucket. *The* 

Schütz, S., Piesche, M., Gorbach, G., Schilling, M., Seyfert, C., Kopf, P., Deuschle, T., Sautter,

Spelter, L. E., Schirner, J. and Nirschl, H., (2011). A novel approach for determining the flow

Spelter, L. E. & Nirschl, H. (2010). Classification of Fine Particles in High-Speed Centrifuges.

Spelter, L. E., Steiwand, A. & Nirschl, H. (2010). Processing of dispersions containing fine

Stahl, S., Spelter, L. E. & Nirschl, H. (2008). Investigations on the Separation Efficiency of Tubular Bowl Centrifuges. *Chemical Engineering & Technology*, 31, 11, 1577-1583

Stevens, A. B. & Hrenya, C. M. (2005). Comparison of soft-sphere models to

Taylor, A. R. (1946). Concentration of the Rabbit Papilloma Virus with the Sharples

Trawinski, H. (1959). Die äquivalente Klärfläche von Zentrifugen. *Chemiker Zeitung* 

van Wachem, B. G. M. & Almstedt, A. E. (2003). Methods for multiphase computational

Wang, B. & Yu, A. B. (2006). Numerical study of particle–fluid flow in hydrocyclones with

Wardle, K. E., Allen, T. R. & Swaney, R. (2006). Computational Fluid Dynamics (CFD) Study

Wilcox, D. C. (1988). Reassessment of the scale-determining equation for advanced

Yakhot, V. & Orszag, S. A. (1986). Renormalization-group analysis of turbulence. *Physical* 

Zubkov, V. A. & Golovko, Y. D. (1968). Liquid Flow in the Extraction Head of a Tubular-Centrifuge Rotor. *Khimicheskoe i Neftyanoe Mashinostroenie*, 12, 18-19

of the Flow in an Annular Centrifugal Contactor. *Separation Science and Technology*,

Supercentrifuge. *Journal of Biological Chemistry*, 163, 1, 283-287

fluid dynamics. *Chemical Engineering Journal*, 96, 1-3, 81-98

turbulence models. *Aiaa Journal*, 26, 11, 1299-1310

different body dimensions. *Minerals Engineering*, 19, 1022-1033

N., Popp, E. & Warth, T. (2007). CFD in der mechanischen Trenntechnik. *Chemie* 

patterns in centrifuges by means of Laser-Doppler-Anemometry. Chemical

particles or biological products in tubular bowl centrifuges. *Chemical Engineering* 

measurements of collision properties during normal impacts. *Powder Technology*,

Centrifuge. *Chemical Engineering & Technology*, 32, 5, 719-725

Sokolow, W. J. (1971). *Moderne Industriezentrifugen*, VEB Verlag Berlin, Leipzig

Supercentrifuge. *Analytical Chemistry*, 20, 4, 379-381

*American Society of Mechanical Engineers*, 118, 173-182

*Chemical Engineering Science*, 46, 8, 2143-2152

Engineering Science, In Press, Corrected Proof

Stahl, W. H. (2004). *Industrie-Zentrifugen*, DrM Press, Männedorf

*Chemische Apperatur*, 83, 18, 606-612

*Review Letters*, 57, 14, 1722-1724

*Chemical Engineering & Technology*, 33, 8, 1276-1282

*Ingenieur Technik*, 79, II, 1777-1796

*Science*, 65, 14, 4173-4181

154, 2-3, 99-109

41, 2225-2244


Gösele, W. (1974). Strömungen in Überlauf-Zentrifugen und ihr Einfluß auf die

Hartl, J. & Ooi, J. Y. (2008). Experiments and simulations of direct shear tests: porosity,

Hirt, C. W. & Nichols, B. D. (1981). Volume of Fluid (VOF) method for the dynamics of free

Horanyi, R. & Nemeth, J. (1971). Theoretical Investigation of Clarification Process in a Tube Centrifuge. *Acta Chimica Academiae Scientarium Hungaricae*, 69, 1, 59-75 Kelecy, F. J. (2008). Coupling momentum and continuity increases robustness. *ANSYS* 

Kynch, G. J. (1952). A theory of sedimentation. *Transactions of the Faraday Society*, 48, 2, 166-

Landman, K. A. & White, L. R. (1994). Solid/liquid separation of flocculated suspensions.

Launder, B. E., Reece, G. J. & Rodi, W. (1975). Progress in development of a Reynolds-Stress

Launder, B. E. & Spalding, D. B. (1974). The numerical computation of turbulent flows.

Li, Y., Zhang, J. P. & Fan, L. S. (1999). Numerical simulation of gas-liquid-solid fluidization

Li, Y. J., Xu, Y. & Thornton, C. (2005). A comparison of discrete element simulations and

Lun, C. K. K. & Savage, S. B. (1987). A simple kinetic theory for granular flow of rough,

Mainza, A., Narasimha, M., Powell, M. S., Holtham, P. N. & Brennan, M. (2006). Study of

Mindlin, R. D. (1949). Compliance of elastic bodies in contact. *Journal of Applied Mechanics-*

Mousavian, S. M. & Najafi, A. F. (2009). Numerical simulations of gas-liquid-solid flows in a

Ni, L. A., Yu, A. B., Lu, G. Q. & Howes, T. (2006). Simulation of the cake formation and

Nowakowski, A. F., Cullivan, J. C., Williams, R. A. & Dyakowski, T. (2004). Application of

Perardi, T. E. & Anderson, N. G. (1970). K-series centrifuges : III. Effect of core taper on

Perardi, T. E., Leffler, R. A. A. & Anderson, N. G. (1969). K-series centrifuges II. Performance

CFD to modelling of the flow in hydrocyclones. Is this a realizable option or still a

hydrocyclone separator. *Archive of Applied Mechanics*, 79, 5, 395-409

growth in cake filtration. *Minerals Engineering*, 19, 10, 1084-1097

particle capture efficiency. *Analytical Biochemistry*, 34, 1, 112-122

research challange? *Minerals Engineering*, 17, 661-669

of the K-II rotor. *Analytical Biochemistry*, 32, 3, 495-511

systems using a combined CFD-VOF-DPM method: bubble wake behavior.

experiments for 'sandpiles' composed of spherical particles. *Powder Technology*, 160,

inelastic, spherical particles. *Journal of Applied Mechanics-Transactions of the Asme*, 54,

flow behaviour in a three-product cyclone using computational fluid dynamics.

turbulence closure. *Journal of Fluid Mechanics*, 68, APR15, 537-566

*Computer Methods in Applied Mechanics and Engineering*, 3, 2, 269-289 Leung, W. (1998). *Industrial Centrifugation Technology*, McGraw.Hill, New York

contact friction and bulk friction. *Granular Matter*, 10, 4, 263-271

boundaries. *Journal of Computational Physics*, 39, 1, 201-225

Sedimentation. *Chemie Ingenieur Technik*, 46, 2, 67

*Advances in Colloid and Interface Science*, 51, 175-246

*Chemical Engineering Science*, 54, 21, 5101-5107

*Minerals Engineering*, 19, 10, 1048-1058

*Transactions of the Asme*, 16, 3, 259-268

*Advantage*, 2, 2, 49-51

176

3, 219-228

1, 47-53

Reuter, H. (1967). Sedimentation in der Überlaufzentrifuge. *Chemie-Ing. Tech.*, 39, 9, 548-553


**CFD and Thermography Techniques** 

**Applied in Cooling Systems Designs** 

The focus of this work consists in optimizing the water flow inside a water cooled electric motor frame, aiming at the maximization of the power/frame size ratio, the minimization of pressure drop and the avoidance of hot spots. For the development of this work

For many years water cooled electric motors have been used by industry, especially in specific applications that require features that can not be provided by conventional fan cooled motors. Among theses features are a better power/frame size relation, low noise

Cooling by means of water circulation is frequently used in large motors, typically frame sizes above IEC 315. A justification for this practice not to be widely used in smaller motors is its higher cost compared to air cooling systems. However, provided that it presents some technical and/or economical viability, water cooling can also be used in smaller frame sizes. The technological progress have resulted in the development of new tools that facilitate the study of this motor type, such as the computational fluid dynamics and the thermography techniques, which consist, respectively, in the use of numeric models applied to fluid

The following advantages and disadvantages of this kind of motor can be mentioned, in

computational fluid dynamics (CFD) and thermography techniques were used.

mechanics and the obtaining of the motor surface temperature distribution.

The outside waste deposit on the frame does not damage the motor cooling;

 The heat removed from the motor is not directly dissipated into the environment; Possibility of using the same cooling fluid of the load and/or machine where it is

No problems with waste deposit on the external fan blade;

It can be used in totally enclosure environment;

**1. Introduction** 

levels, enclose applications, among others.

**2. Benefits and drawbacks** 

**2.1 Advantages** 

installed;

 Lower noise level; Higher efficiency;

comparison with air cooled motors.

Greater power/frame size ratio;

Samuel Santos Borges and Cassiano Antunes Cezario *Research and Technological Innovation Department, WEG* 

*Brazil* 

Zubkov, V. A. & Golovko, Y. D. (1969). Flow of Liquid in Rotor Outlet of a Tubular Centrifuge. *International Chemical Engineering*, 9, 3, 403-406 **7** 
