Computer simulations of sodium formate solution in a mixing tank

Computer simulations of sodium formate solution in a mixing tank

Communications in Nonlinear Science and Numerical Simulation 13 (2008) 380–399 Computer simulations of sodium formate s...

2MB Sizes 0 Downloads 3 Views

Communications in Nonlinear Science and Numerical Simulation 13 (2008) 380–399

Computer simulations of sodium formate solution in a mixing tank S.M. Mousavi


, P. Zamankhan b, A. Jafari




Department of Chemical and Petroleum Engineering, Sharif University of Technology, Tehran, Iran Department of Energy and Environmental Engineering, Lappeenranta University of Technology, Lappeenranta, Finland Received 23 February 2006; received in revised form 1 May 2006; accepted 1 May 2006 Available online 22 June 2006

Abstract Traditionally, solid–liquid mixing has always been regarded as an empirical technology with many aspects of mixing, dispersing and contacting were related to power draw. One important application of solid–liquid mixing is the preparation of brine from sodium formate. This material has been widely used as a drilling and completion fluid in challenging environments such as the Barents Sea. In this paper, large-eddy simulations of a turbulent flow in a solid–liquid baffled cylindrical mixing vessel with large number of solid particles are performed to obtain insight into the fundamental aspects of a mixing tank. The impeller-induced flow at the blade tip radius is modeled by using the dynamic-mesh Lagrangian method. The simulations are four-way coupled, which implies that both solid–liquid and solid–solid interactions are taken into account. By employing a soft particle approach the normal and tangential forces are calculated acting on a particle due to viscoelastic contacts with other neighboring particles. The results show that the granulated form of sodium formate may provide a mixture that allows faster and easier preparation of formate brine in a mixing tank. In addition it is found that exceeding a critical size for grains phenomena, such as caking, can be prevented. The obtained numerical results suggest that by choosing appropriate parameters a mixture can be produced that remains free-flowing no matter how long it is stored before use.  2006 Elsevier B.V. All rights reserved. PACS: 47.11.j; 47.27.Ep; 47.51.+a Keywords: Large-eddy simulation; Sodium formate; Mixing tank; Dissolution; Fluid–structure interaction; Particle dynamic model

1. Introduction Stirred tanks have wide applications in the chemical and process industries. Crystallization operations [1], liquid–liquid extractions [2], biological fermentations [3], and heterogeneous catalytic reactions [4] are just a few examples of industrial operations usually carried out in mixing tanks. In the aforementioned applications, * Corresponding author. Address: Department of Chemical and Petroleum Engineering, Sharif University of Technology, Tehran, Iran. Tel.: +98 21 66165494; fax: +98 21 66005417. E-mail address: [email protected]edu (S.M. Mousavi).

1007-5704/$ - see front matter  2006 Elsevier B.V. All rights reserved. doi:10.1016/j.cnsns.2006.05.002

S.M. Mousavi et al. / Communications in Nonlinear Science and Numerical Simulation 13 (2008) 380–399


mixing tanks are required to fulfill several needs including blending of miscible liquids, dispersion of gases or immiscible liquids into a liquid phase, suspension of solid particles, heat and mass transfer enhancements, and chemical reactions. Among the various industrial unit operations involved with multi-phase systems, agitation of solid–liquid systems are commonly encountered in catalytic reactions, leaching, and polymerization. Despite their widespread use, the complex 3-D recirculation and turbulent flows in the tank makes designing and optimizing of the reactor usually limited to pilot plant tests and empirical formulation, even for single phase applications. A number of investigations have been published on the fluid dynamic properties of solid liquid systems in mixing tanks for achieving empirical information. Most investigations, which are focused, on the distribution of solid particles are given in Baldi et al. [5], Barris and Baldi [6], MacTaggart et al. [7], and Nasr-El-Din et al. [8]. In addition, attempts with the criteria of suspension include Chudacek [9], and Zwietering [10]. In spite of the importance and variety of applications, the problem of designing and scaling up mixing tanks has been tackled mainly by means of semi-empirical methods. As suggested by Tatterson [11], economical losses are huge as a consequence of uncertainties in the design of mixers. It is believed that a significant improvement on mixing tank design can be obtained by developing numerical models, which take into account the real flow field inside the vessel. Efforts have been made during the last years towards the development of predictive methods based on computational fluid dynamics (CFD) and capable of providing detailed information on flow and turbulence fields of mixing tanks [12]. Sufficiently strong rotation in mixing tanks is observed to considerably change the turbulence dynamics leading to highly anisotropic structuring of turbulent eddies [13]. On the other hand, CFD is incapable of solving all flow scales in industrial applications, within a reasonable amount of time and memory space. A direct numerical simulation (DNS) of the turbulent flow at industrially relevant Reynolds numbers (Re P 104) would require an enormous amount of grid cells and time steps to capture all relevant time and length scales in the flow. In a large-eddy simulation (LES), the smallest scales are modeled with a subgrid-scale model, thereby reducing the amount of computational resources [14]. In this light, the advanced modeling techniques such as Large-eddy Simulation (LES) would be required in order to predict the flow dynamics in a mixing tank correctly. Recent applications of LES to the hydrodynamic simulation of industrial equipment such as baffled stirred tanks are presented in [15–19]. Large-eddy simulation (LES) is a simulation technique [14] based on the decomposition of the flow field into large and small-scale structures. In this case the large-scale structures are directly simulated in three-dimensional, time dependent mode, whereas the feed-back effects of the small-scale structures are modeled. Preparation of brine from sodium formate is one of the important applications of solid–liquid mixing. The sodium formate, as illustrated in Fig. 2, consists of white hygroscopic crystalline odorless granules,

Fig. 1. Schematic and some dimensions of standard baffled tank with Rushton radial turbine used in the present attempt. (a) Gas–liquid mixing tank in which the interface is located at H = 260 mm. (b) Three phase mixing tank with the solid phase is shown to be deposited at the bottom. (c) Different views of the 6 impellers Rushton turbine with its characteristic dimensions.


S.M. Mousavi et al. / Communications in Nonlinear Science and Numerical Simulation 13 (2008) 380–399

Fig. 2. (a) Sodium formate described as white granules with molecular weight of 68.01, and chemical formula of HCOONA. (b) A sample of sodium formate used in the computer simulation reported in this study.

which has been widely used to prepare drilling and completion fluids. It is often used in challenging environments such as the Barents Sea. Formate brine solutions including those made from sodium are used in the oilfield as drilling fluids where fluid viscosity must be low. These fluids are made up of the formate solution itself and water soluble polymers. The fluids are colorless and can be used to drill into the pay zone without threat of plugging or damaging formations. Sodium formate chemically reduces other components by donating an electron or electrons. Formic acid and oxalic acid are prepared from this material. Sodium formate is used in the manufacture of sodium hydrosulfite, a common reductive bleaching chemical. More applications of sodium formate such as that in the recovery of precious metals from acidic effluents have been discussed in Julsing and McCrindle [20]. Sodium formate is also used to improve the brightness and color in printing fabrics and paper. The process using sodium formate powder is a difficult one because the individual particles of powder can stick together during storage and fuse into a solid mass that is difficult to break up and hydrate in a brine-mixing tank. Published data for solubilities of sodium formate in water include Groscuff [21] and Sidgwick and Gentle [22]. In the light of above discussion, in the present attempt an overview is presented of collective processes in liquid–particle flows useful for developing a simplified model for particle dynamic type simulations of dense liquid–particle flows. The governing equations of the liquid phase with air above are solved using LES technique. The particle motion is predicted by a Lagrangian method. Particles are assumed to behave as visco-elastic solids during interactions with their neighboring particles and inter-particle normal and tangential contact forces between particles are calculated using a generalized Hertzian model [23]. The other forces on a particle that are taken into account are gravitational, and drag force resulting from velocity difference with the surrounding liquid. A simulation of gas–particle flow is performed for predicting the flow dynamics of dense mixture of liquid and particles in a mixing tank.

S.M. Mousavi et al. / Communications in Nonlinear Science and Numerical Simulation 13 (2008) 380–399


2. Theoretical 2.1. The LES model The numerical modeling of these complicated mixing processes is a daunting task. Direct numerical simulation (DNS) provides the most exact approach in which the mechanism involved in turbulent mixing can be accurately represented. DNS requires resolving the smallest eddies which makes the approach prohibitively expensive even with the most powerful computers of the present day, and foreseeable future as well. On the other hand, the popular approaches based on the Reynolds-averaged Navier–Stokes (RANS) equations amount to averaging out the large eddies that are mainly responsible for mixing. Recently, the large-eddy simulation (LES) approach has been used as an intermediate method between the extremes of DNS and RANS. The premise underlying LES is that any turbulent flow consists of eddies, and involves a wide spectrum of scale. In LES, large eddies are resolved directly, and small eddies are modeled. Large eddies are more dependent on the geometries and boundary conditions of the flow involved, and are highly anisotropic. Small eddies are less dependent on the geometry, and tend to be more isotropic, and consequently are more universal. The chance of finding an universal model is much higher when only small eddies are modeled. In the LES approach, the governing equations are obtained by spatially filtering the Navier– Stokes equations. The large turbulent scales are computed explicitly, while the small scales are modeled using one of a number of available subgrid-scale (SGS) models. The SGS models describe interactions between the resolved and unresolved scales. The LES approach is more general than the RANS approach, and avoids the RANS dependence on boundary conditions for the large-scale eddies. Like DNS, the LES approach gives a three-dimensional, time dependent solution. The required computational resources for LES are between those of the DNS and RANS approaches. The LES model can be used at much higher Reynolds numbers than DNS because the computational effort is independent of the Reynolds number if the small scales obey the inertial range spectrum and the near wall effects are not important. 2.2. Governing equations Turbulence in a steadily rotating fluid, such as found in brine-mixing tanks, plays an important role in solid dispersion. As illustrated in Fig. 3, sufficiently strong rotation is observed to considerably alter the turbulence dynamics, in particular leading to anisotropic structuring of turbulent eddies. The importance of rotation for the large scales of turbulence is determined by the reciprocal of the turbulent Rossby number, Ro, where it is the ratio of a rotation time scale, X1, to a characteristic time, L/u 0 , for evolution of the large scales of the turbulence in the absence of rotation. Hence, a large Rossby number represents a case with small effects of rotation and if it is sufficiently high one may neglect rotation altogether. On the other hand, a small enough Rossby number implies strong rotation and, therefore, not enough time for the large scales of turbulence to act on themselves significantly during a rotation period. 2.2.1. Air–sodium formate solution flows A rotating air–sodium formate solution (which is called solution after that) system as illustrated in Fig. 1(a) may be described using the volume-of-fluid (VOF) technique [24]. This technique provides an excellent approximation when the ratio of liquid to gas densities is large. In the present study, the VOF is used for modeling flows of air and solution as immiscible fluids with interfaces effects in a brine-mixing tank. Here, the computations are performed on a Cartesian grid, with the interface (between solution and air, in which air can only apply a pressure on solution) being localized by calculating the fraction of each computational cell occupied by one of the two phases. Surface tension effects may be incorporated by modeling the capillary stresses as equivalent body force acting in the immediate vicinity of the interface. In this case, an algorithm is required to track the surface as a sharp interface moving through a computational grid. In the VOF, the velocity of air and solution is considered equal at the interface. Recall that the LES has been successfully applied for simulations of the flow driven by a Rushton turbine [15–19]. Hence, in the current attempt the aim is to combine the best features of the VOF method with those of the LES in order to achieve more accurate simulations of air–solution flows in a brine-mixing tank.


S.M. Mousavi et al. / Communications in Nonlinear Science and Numerical Simulation 13 (2008) 380–399

Fig. 3. (a) The instantaneous velocity vector field in a cutting xz-plane passing through the rotor axis of a mixing tank, whose schematic is illustrated in Fig. 1(a), at the tank Reynolds number (defined as ND2/m) of Retank = 48500. (b) The interface of air and water in the mixing tank is shown to obtain better visualization. (c) A high resolution version of the velocity vector field in a cutting xy-plane lactated below of the impellers.

The filtered continuity, momentum for an isothermal air–solution system may be given as o ð ui Þ ¼ 0: oxi o o o p o rij osij ðqm  ðqm  þ þ  qm gdi3 þ F sv ui Þ þ ui  uj Þ ¼  i ; ot oxi oxi oxj oxj and

ð1Þ ð2Þ

  2 oui ouj rij ¼  lm vm;m dij þ lm þ ; 3 oxj oxi

where ui is liquid velocity, qm is mixture density, p is pressure, rij is viscous stress tensor [25], g is acceleration due to the gravity, F sv i is volume form of the surface tension, lm is mixture viscosity, and dij is Kronecker delta. Overbar in this study represents spatial filter. In the present study, the unresolved part of F sv i is neglected.  2 Þ1=2 expðcjxi  ni j2 =D  2 Þ, the subgrid-scale Applying a Gaussian spatial filter such as Gðxi  ni Þ ¼ ðc=pD stress tensor may be modeled as sij ¼ C ij þ Rij ¼ ui uj  ui uj . Here, C ij ¼ ui u0j þ uj u0i represents the interaction between large and small scales, and Rij ¼ u0i u0j represents the interaction between the subgrid scales. Note that the  ui  uj term in Eq. (2) requires a second application of the filter. To remedy this further decomposition is required  ui  ui  ui  ui  uj ¼ ð ui  uj   uj Þ þ  uj ¼ Lij þ  uj ;


where Lij indicates interactions among the large scales. Using Eq. (3), the filtered momentum equation may be given as   ^ o ui o o p o o ui o uj o s ij þ ð ui  þm þ ; ð4Þ uj Þ ¼   oxi oxj oxj oxi ot oxj oxj ^

where s ij ¼ Lij þ C ij þ Rij ¼ ui uj   ui  uj , and m is kinematic viscosity of the liquid.

S.M. Mousavi et al. / Communications in Nonlinear Science and Numerical Simulation 13 (2008) 380–399


In order to model the subgrid-scale tensor, Menon and Kim [26], suggested that based on the similarity assumption of different scale sij = CkLij, a kinetic energy ktest can be defined of a test level that includes all b given as energy between two length scales ðD < l < DÞ, 1 b b k test ¼ ð ud ð5Þ k uk  uk uk Þ: 2 Assuming a similar representation of these two levels, expressions may be obtained for subgrid stress tensor 1=2

sij ¼ 2C s Dk SGS S ij þ dij =3skk ;

ð6Þ b 1=2 Scij .  Dk test

where C s ¼ 1=2Lij rij =rpq rpq and rij ¼ Dissipation of energy at the test level takes place between the two length scales may be expressed as ! oud o ubi o ubi k 3=2 i oui etest ¼ ðm þ mt Þ  ¼ C e test : ð7Þ b oxj oxj oxj oxj D Thus, eSGS ¼ C e

3=2 k test : D

The model of Menon and Kim [26], gives a transport equation for kSGS that uses the local value of Cs in the diffusion and Ce in the dissipation term. That is   3=2 ok SGS oui k SGS o ok SGS k 1=2 þ ¼ ðC s Dk SGS þ mÞ ð8Þ þ 2mt S ij S ij  C e SGS ; oxi ot oxi oxi D 1=2

where mt ¼ C s Dk SGS . Treating gas in the tank as incompressible greatly simplifies the analysis. However, a careful treatment of free surface is required to avoid producing any incorrect motion of the surface since it is assumed to move with the average velocity of air and solution. Following Chen et al. [27], the density and viscosity in the momentum Eq. (2) may be expressed as qm ¼ qg ð1  F Þ þ F qsol ; lm ¼ lg ð1  F Þ þ F lsol ;


where F (0 6 F 6 1) is fraction of a control volume occupied by the liquid, qg is gas density, qsol is the solution density, lg is gas viscosity, and lsol is solution viscosity. Combining the interface transfer equation in terms of F, namely oF/ot + viF,i = 0, with continuity (1), gives oF þ ðvi F Þ;i ¼ Fvi;i : ot


Here, F is used to track the interface instead of a rapid change in density. The unit normal vector for the interface is given by ni = F,i/jF,ij. This scheme results in somewhat sharp interfaces localized in one control volume. However, the interface having surface tension may not be necessarily aligned with logical mesh coordinates. In this light, the interface in the computational domain may only be represented as a finite thickness transition region comparable to mesh spacing through which the liquid volume fraction varies from zero to one. An expression for the volume form of the surface tension force, F sv i , is given by [25] F ;i ; ð11Þ F sv i ¼ rj ½F  r is surface tension coefficient and j represents the curvature of the interface and is defined as   1 F ;j jF ;i j;j  F ;ii : j¼ jF ;i j jF ;i j where F sv i is non-zero only within the transition regions where the liquid volume fraction varies smoothly from zero to one. The square brackets in Eq. (11) represent the jump in the value of F across the interface.


S.M. Mousavi et al. / Communications in Nonlinear Science and Numerical Simulation 13 (2008) 380–399

The large-eddy simulations discussed in this section are performed for a mixing tank, which contains solution with air above, with a standard Rushton turbine, as illustrated schematically in Fig. 1(a). The focus is limited to a case with the Reynolds number of 48,500. In the present attempt a sliding mesh method is used by which flow pattern can be calculated without the use of any experimental boundary conditions. The tank is divided to two zones, which are the impeller zone and the tank zone. The latter includes the tank wall, the tank base and the baffles. The unstructured grid in the impeller zone rotates with the impeller, whereas the grid in the tank remains stationary. The aforementioned grids slide past each other at a cylindrical sliding interface. No-slip boundary conditions are employed at the impeller blades, the baffles, and the tank walls. The governing equations listed in this section are fully solved based on control volume technique. A similar method as that detailed in [27] is implemented for the interface tracking. The sample results as presented in Fig. 3 are obtained using a total number of grid nodes that is more than 106. The central differencing scheme is used for spatial discretization and the time is advanced via a second-order implicit scheme. 2.2.2. Particle–solution flows Solid particles. In general, solid particles such as those shown in Fig. 2(b) with mass m in a liquid obey the Langevin equation given as m

dV pi ¼ F fi þ F gi þ F ii þ F bi ; dt


where V pi is particle velocity, F fi is drag force from fluid, F gi is gravitational force, F ii is force due to particle– particle interaction, and F bi is Brownian force due to thermal motion of the water. To simplify the equation of motion for the solid particles it would be useful to define the dimensionless quantities given as below d p jui  V pi j ; mw !1 2 qw d p jui  V pi j St ¼ ð1  /s Þ ; 9 mw qp

Rep ¼ ð1  /s Þ



Fr ¼

ðð1  /s Þjui  V pi jÞ ; gL

where Rep is particle Reynolds number, /s is solid volume fraction of sodium formate (solid phase), dp is particle diameter, St is Stokes number, qw is water density, qp is material density of the grains, and L is linear size of the tank. Fr, Froude number, is the ratio of the total kinetic energy of the particle to the gravitational potential. Using the data presented in Table 1, the order of magnitude of dimensionless quantities listed in (13) for the mixing tank as illustrated in Fig. 1, are estimated as Rep  100, St  100, and Fr  1. Since the work done by the drag force over the size of particle is much larger than thermal energy, the random force F bi can be neglected. Thus for the case Rep  100, and St  100 with nonlinear, visco-elastic, particle–particle interactions, the Langevin equation for translational motion of the jth particle in the mixing tank may be reduced to j X dV ji F D ðnÞ ¼ i  gdiz þ F i;jp ; dt mp p¼1



where F D i is drag force, mp represents particle mass, and Nj is number of neighboring particles in contact with the jth particle at time t. An expression for the drag force on particles with both the translational and rotational motion, as illustrated in Fig. 4, may be given by [28]   1 xR p p p D F i ¼ qw jui  V i jA C D ðui  V i þ C LR ðui  V i Þ  ð15Þ þ F LG ; 2 jxR j where A is projected area of particle, CD is drag coefficient, CLR is lift coefficient due to particle rotation, xR is particle rotational velocity relative to liquid vorticity, and FLG is lift force due to velocity gradient.

S.M. Mousavi et al. / Communications in Nonlinear Science and Numerical Simulation 13 (2008) 380–399


Table 1 Physical and chemical properties of materials used in the simulation Material


Value Reported

Sodium formate

Used in the model

1.258 at 20 C (solution) 1.92 at 20 C (solid phase) 97 g/100 g water at 20 C White, deliquescent, granules Spherical particles as or crystalline powder illustrated in Fig. 2(b) Not reported 0.65 (measured) Not reported 0.244 (assumed) Not reported 2.53 · 1010 kg/m s (assumed) Not reported 0 kg/m s (assumed) Not reported 9.87 · 106 s (assumed) Not reported 6.3 · 1010 Pa (assumed) HCOONa 68.01

Relative density Solubility Material type Surface friction coefficient Poisson’s ratio Instantaneous shear modulus Long time shear modulus The relaxation time Elastic modulus Chemical formula Molecular weight


Temperature Kinematic viscosity Density

20 C 106 m2/s 998.2 kg/m3


Temperature Kinematic viscosity Density

20 C 1.5 · 105 m2/s 1.2 kg/m3


Poisson’s ratio Elastic modulus Density

0.29 2.1 · 1011 Pa 7890 kg/m3

The sodium formate particles have surface roughness on many different length scales. When two particles are brought into contact, the area of contact will only be a fraction of the nominal contact area. The contact regions are small areas where asperities from one ball are squeezed against those of the other ball. By assuming that the asperities will deform elastically with a finite viscous relaxation time [29], s, an expression may be obtained for the contact force per unit mass acting on the jth particle due to contacts between this particle and its neighboring pth particles at time t, given as [23] _ ! _ _ dd K s jp n 3=2 3=2 1=2 2 F i;jp ¼ ð2E=ð3ð1  m2p Þmp ÞÞd p d jp þ ðG0 =ðG0  G1 ÞÞd jp þ ct V imp ð16Þ k i;jp  ðk t vjp l;jp el Þei ; mp dt _

where E is Young modulus, mp is Poisson’s ratio, d jp is deformation of largest asperities (which represents overlapping between the pth and jth particles), Kn is coefficient characterizing the viscous behavior of the grains as given in Table 1, G0 is instantaneous (glassy) shear modulus, G1 is long time shear modulus, kt is elastic modulus in a tangential direction, Following Silbert et al. [30] the tangential displacement is denoted by vjp, and el is unit vector. The surface asperities are responsible for the tangential forces acting between the colliding pairs. The first and second terms on right hand side of (16) represent the normal and tangential components of the contact V imp V imp force, respectively. In Eq. (16) ei is the tangential unit vector defines as, ei ¼ eimq k m;jp ðeqst es;jp k t;jp Þ, where, es;jp , V imp imp ¼ V imp =jV j. In this case, represents the unit vector in the direction of relative impact velocity defined as es;jp s;jp s;jp Coulomb friction law describes friction between two colliding grains with a surface friction coefficient, lp, when there is mutual slipping at the point of contact. The magnitude of vjp, which represents the tangential displacement, is calculated as necessary to satisfy jFtjpj = lpjFi,njpki,jpj. Otherwise, the contact surfaces are considered as stuck while jFtjpj < lpjFnjpj. jFtjpj and jFnjpj are the magnitude of the normal and tangential components of the contact force, respectively. Here, the rate of change of tangential displacement, vjp, is given by dvjp =dt ¼ V imp i;jp ei . The displacement, vjp, should be set initially to zero when a new contact is established and once the contact is broken, all memory of the prior displacement will be lost. As illustrated in Fig. 5, frictional forces induce torques on particles. The


S.M. Mousavi et al. / Communications in Nonlinear Science and Numerical Simulation 13 (2008) 380–399

Fig. 4. The complex velocity vector field around rotating grains of solid particles in a solution. In this case, the drag force includes the lift force due to particle rotational velocity relative to liquid vorticity, which is different than the lift force due to velocity gradient. Vectors are color coded by their velocity magnitudes, where the red is for highest velocity and blue represents the lowest. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

PN torque on particle j may be defined as [30] T i;j ¼  12 p¼1 eiqr mj rjp k q;jp F r;jp . In some cases, the rotation rates of turbine, N, as large as 2 Hz is induced due to frictional contact between rough grains of sodium formate. Hence, Eq. (14) must be augmented by a torque equation for the rotational motion of particle j which is given as Jj

dxi ¼ T i;j : dt


where Jj is moment of inertia of particle j, xi is rotational velocity of particles, and Ti,j is torque on particle j. Assuming that the apparent surface of contact is built up of a lager number of hierarchically ordered asperities, Brilliantov et al. [31] suggested that the friction coefficient, lp, may be expressed in terms of mesoscopic parameters. However, in the present attempt the friction coefficient is found by comparing the geometry of a dry spill of sodium formate particles observed experimentally and that obtained using simulations. A detailed description of the method for estimating the coefficient of friction is presented in the following. To obtain an expression appropriate to spherical particle-flat wall contact from Eq. (16) r has to replace 2r.

S.M. Mousavi et al. / Communications in Nonlinear Science and Numerical Simulation 13 (2008) 380–399


Fig. 5. Temporal evolution of the velocity of colliding particles with diameter of dp = 5 mm. The particles are color-coded using their velocity magnitudes. Instantaneous configurations in (a) and (b) are each separated by t  106 s. Coefficient of friction. This section addresses spills of sodium formats onto a flat surface which collect in a pile. Molecular dynamic type simulation techniques may provide insights into dynamics of the pile, which are not very well understood. Traditional molecular dynamics analysis relies on the classical Hertzian impact theory and uses particle contacts of finite duration. In this case, a model for collision between particles was used by artificially increasing the duration of a contact between the particles, as suggested by Silbert et al. [30]. However, the disadvantage of the traditional approach is that the softer, normal interactions between the particles have to be assumed. Hence, prediction of the contact network composed of frictional contacts by the traditional model should be considered with caution. By extending the existing formalism an accurate contact dynamics algorithm is proposed by which the collective intermittent dynamics found in piles may be captured. Here, the results of a molecular dynamic type simulation of spills of rough sodium formate whose physical properties are listed in Table 1, consisting of 507 mono-sized spherical particles, are presented using the generalized model discussed in the preceding section. The resulting shape of spill of glass beads may be controlled by the distance between the leak and the surface where the leak collects, and the smoothness of the surface on which the material collects. In this light, the simulation includes the effect of static friction, which is found to be crucial in maintaining a stable spill. Moreover, the diameter of beads are set to r = 1 mm, which is large enough that no aggregates can be produced via Van der Waals forces. As illustrated in Fig. 6(a), which represents the result of experimentation using rough glass beads with diameter r = 1 mm, a typical spill exhibits different distinct parts including a rounded shape on top, a linear region characterized by the angle of repose and a tail at the bottom. The angle repose, which is the angle made between the surface and the outside of the cone, achieved in a system as shown in Fig. 6(a) seems to be slightly less than the maximum angle of repose produced using simulation techniques detailed in the preceding section.


S.M. Mousavi et al. / Communications in Nonlinear Science and Numerical Simulation 13 (2008) 380–399

Fig. 6. (a) A snapshot of a spill of roughly mono-sized sodium formate granules onto a flat surface, which collect in a right circular conical shape. The particle diameter is dp  1 mm. The straight line represents the slope of the linear region. In this case, the angle of repose is roughly 33. (b) An instantaneous configuration of particles slightly before the 507th particle collides in the pile.

Note that the ease with which spilled material can move along the collecting surface, minor grooves, ridges, and hills and valleys on a collecting surface may slightly affect the spill geometry. Fig. 6(b) represents the resulting shape of spill of glass beads using simulations where a static coefficient of friction, lp = 0.65, is used. In order to obtain an angle of repose of hf  33, the coefficient of static friction is set to 0.65 in the present simulations. A slightly higher coefficient of static friction has to be used to model the sliding movements of particles over a flat surface made of steel. Here, the coefficient of static friction is set to 0.8 for grain steel sliding contacts. The base of the flat surface is supported so that no movement is possible in any direction and the effects of the gas are neglected. The experimental results are fairly well reproduced by the present model implying the correct physical framework has been incorporated. The linear region characterized by the angle of repose as depicted in Fig. 6(a) is satisfactorily reproduced by model detailed in the preceding section whose results are shown in Fig. 6(b). Solubility. Solubility is the amount of sodium formate that will dissolve in the water. Solubility of sodium formate has special significance in design of brine-mixing tanks. In this section, a classical particle dissolution rate expression is developed, which will be used to predict particle dissolution rate phenomena in the discrete model detailed in preceding sections. Consider a single grain of sodium formate, which dissolves in the surrounding water, as shown in Fig. 7. When there is little sodium formate already in solution, dissolving takes place rapidly. As the solution approaches the point where no sodium formate can be dissolved, dissolving takes place more slowly. In this case, adding heat is one way to increase solubility. The goal is to find both the dissolution rate and the concentration profile around a single, stationary grain of sodium formate. The mass balance in spherical coordinates originating from the center of grain may be given as

S.M. Mousavi et al. / Communications in Nonlinear Science and Numerical Simulation 13 (2008) 380–399


Fig. 7. A diffusion model for the dissolution of a single grain of sodium formate in water. The concentration field within the solid phase Cs(r) = Cs = constant (r 6 Rp). The concentration jumps across the interface, DCeq = Cl0  Cs, at the solid surface. The solution has lower average concentration of sodium formate, C1, than the equilibrium value Cl0.

    oC l 1 o 2 oC l Dsol r ¼ 2 r or ot or


where r represents the distance from center of particle, and Cl is concentration of sodium formate within the boundary layer defined as qlSF/SF + qw(1  /SF). /SF is the volume fraction of sodium formate (liquid phase) around a grain, and qlSF is the sodium formate density. Assuming that the coefficient of diffusivity of sodium formate (liquid) in water, Dsol, changes insignificantly in the radial direction, Eq. (18) reduces to  2   oC l o C l 2 oC l ¼ Dsol þ ð19Þ ðr > Rp Þ: r or ot or2 The solution to this problem, using a spherical transform variable u = Clr and then taking Laplace transform with respect to time, may be give as ! Rp ðr=Rp Þ  1 C ¼ erfc pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð20Þ ðr=Rp P 1Þ: r 4Dsol t=Rp where C* is dimensionless concentration, and Rp represents particle radius. The solution (20) indicates that at short times, namely t  R2p =pDsol , an expression for the interfacial flux may be given as J ðtÞ ffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffi Dsol =ptðC l0  C 1 Þ, whereas at long times, namely t  R2p =pDsol , it simplifies to a value, which depends on the grain size, as given by J ðRp Þ ffi Dsol ðC l0  C 1 Þ=Rp . This suggests that using the smaller size grains may increase the interfacial flux. Here Cl0 is concentration of sodium formate at the solid surface, and C1 shows average concentration of sodium formate defined as qlSF/SF1 + qw(1  /SF1), here /SF1 is the volume fraction of sodium formate in mixing tank. The interface as illustrated in Fig. 7 moves during the dissolution p processes. The position of the interface ffiffiffiffiffiffiffiffiffiffiffiffi can be tracked by the interface characteristic K sphere ¼ ðRp ðtÞ  Rp0 Þ= 4Dsol t. It is straight forward to obtain the expressions for the characteristic equation pffiffiffi 2 2 2ðK sphere Þ ½1  pK sphere eðK sphere Þ erfcðK sphere Þ ¼ S; ð21Þ where S is dimensionless saturation defined as facial speed given as pffiffiffiffiffiffiffiffi dRp j ¼ ðK sphere Dsol Þt1=2 ¼ pffi : dt t

ðC l C 1 Þ . ðC l C S Þ

CS is concentration within the solid phase. The interð22Þ


S.M. Mousavi et al. / Communications in Nonlinear Science and Numerical Simulation 13 (2008) 380–399

The diffusion layer thickness may be predicted using expression as given by dBL S ¼ : K sphere 2ðK sphere Þ2


where dBL is diffusion boundary layer. As shown in Fig. 8, the thickness of boundary layer decreases with increasing the magnitude of saturation, jSj. The ratio between particle size and diffusion layer thickness is an important factor in controlling the shape of particle dissolution profiles. The solution concentration C1 increases with time due to continuous dissolution of sodium formate granules in water. However, the local concentration could be nearly uniform everywhere within the mixing tank due to perfect stirring induced by the Rushton radial turbine. In this light, if the kth grain dissolves atPthe rate of 4pJ k ðtÞR2p , then the total rate of Np accumulation of sodium formate in the tank per unit volume is k¼1 4pJ k ðtÞR2p =V water where Vwater is volume of waterR P in tank. Thus, the average local concentration, C1(t), may be given as C 1 ðtÞ ¼ t Np 4p=V water 0 k¼1 J k ðt0 ÞR2k dt0 . In actual situations as illustrated in Figs. 9 and 10, the dissolution can be complicated by the translational and rotational motion of the particles as well particle–particle interactions. In this case, the mass transfer coefficient may be estimated using Sh ¼ 2kRp =Dsol ¼ ½2ð1  /s Þ þ AReap Scb  [32] where Sc is Schmidt number. Currently, no information exists for the coefficient of A and power a. However, it is reasonable to assume b  0.33 [33]. Here k characterizes the stirring induced by the movements of the grain, which brings fresh portions of the water in contact with the sodium formate, thereby increasing the rate of solution. The simplified diffusion layer model presented in this section provides an approximation for estimating the dissolution rate and the concentration profile around the grain of sodium formate in a brine-mixing tank. Liquid phase. To predict the solution flow field in the mixing tank, a generalized form of the Navier– Stokes equations for the solution interacting with sodium formate granules is used. That is ðqsol ð1  /s Þui Þ;i  ðqsol ð1  /s ÞÞ;t ¼ 0;


ðqsol ð1  /s Þui Þ;t þ ðqsol ð1  /s Þui uj Þ;j ¼ ðð1  /s ÞpÞ;i þ ðð1  /s Þrij Þ;j  ð1  /s Þqsol gdiz þ fi :


The last term on the right hand side of Eq. (25) represents the particle effect on the solution, which is given the Pas Nc sum of all hydrodynamic forces on the particles in a computational cell, namely fi ¼ 1=ðV c ð1  /s ÞÞ s¼1 FD i where Vc is volume of a computational cell and Nc represents number of particles in the cell.

Fig. 8. Variations of dimensionless diffusion layer thickness with saturation. Inset: The typical shape of particle dissolution profile. In the present attempt, the grains are assumed to be dissolving in a partially saturated solution with concentration C1.

S.M. Mousavi et al. / Communications in Nonlinear Science and Numerical Simulation 13 (2008) 380–399


Fig. 9. (a) A top view of initial configuration of sodium formate grains in the tank. (b) The configuration of particles after t  1 s of simulation.

The governing equations of the liquid phase are solved using LES technique. The filtered equations are similar to Eqs. (24) and (25) with the difference that on the left-hand side of Eq. (25) the gradient of the subgridscale stress tensor must be included given by sij,j with sij ¼ qui uj  qui quj =q, where q = qsol(1  /s) is apparent density. In order to achieve closure of the LES equations of the gas phase the subgrid-model is presented in Appendix A. Subgrid-scale model (A.8) was originally proposed for compressible turbulent flow. In the present study, the model (A.8) is tested to clarify its relevance for dense water-sodium formate flows in the mixing tank, for which the fluctuating fields might not be small compared with those of the mean fields. 2.3. Numerical method A sliding mesh method has been used by which flow pattern can be calculated without the use of any experimental boundary conditions. The tank is divided to two zones, which are the impeller zone and the tank zone. The latter includes the tank wall, the tank base and the baffles. The unstructured grid in the impeller zone rotates with the impeller, whereas the grid in the tank remains stationary. The aforementioned grids slide past


S.M. Mousavi et al. / Communications in Nonlinear Science and Numerical Simulation 13 (2008) 380–399

Fig. 10. Different views of a typical configuration of sodium formate granules around the Ruston turbine which rotates at X = 400 rpm. Here only about 104 particles are shown.

each other at a cylindrical sliding interface. No-slip boundary conditions are employed at the impeller blades, the baffles, and the tank walls. To divide the geometry into discrete control volumes, more than 106 nodes and 8 · 105, 3D tetrahedral computational cells were used. In addition, roughly more than 7 · 104 wall triangular elements were used. The convergence was considered to be achieved when the conservation equations of mass and momentum were satisfied, which was considered to have occurred when the normalized residuals became smaller than 5 · 105. The normalization factors used for the mass and momentum were the maximum residual values after the first few iterations. It is worth noting that refinement of the grids did not produce any significant differences in the results. The flow is resolved using 5 · 105 tetrahedral cells. The mesh-spacing of this grid should be larger than the particle diameter. The solution of Eqs. (24) and (25) requires specification of the solids fraction at the appropriate grid nodes. The value of solids fraction is obtained from particle dynamic simulations by counting the volume of particles within each computational cell divided by the volume of the cell. The governing equations are fully solved with a second-order finite volume method on a staggered grid where the time is advanced via a second-order implicit scheme. The time step for the liquid phase equals 105 s. The grain size shrinks with time due to dissolution of sodium formate in the water. The new particle size is calculated every 5 · 104 s. The solution of the gas field requires the computational domain, to be divided into cells. The calculation of the gas is based on the numerical solution of the set of partial differential equations (24), (25) and (A.8).

S.M. Mousavi et al. / Communications in Nonlinear Science and Numerical Simulation 13 (2008) 380–399


3. Numerical results By performing the simulations based on the discrete model outlined in the preceding sections (which contains only a few empirical parameters), an increased understanding of the dynamics of dense liquid–particle flows may be obtained in a mixing tank. Here, the main goal is to address problems such as caking in a brine-mixing tank. The system represents a virtual type of brine-mixing tank. The simulation is performed using pure water with air above in a cylindrical mixing vessel. A Rushton turbine is located in the center of the vessel as shown in Fig. 1(a). The rotor and impellers are made of steel whose physical properties as an elastic material are presented in Table 1. The flow is driven by the rotational motion of the Rushton turbine in a xy  plane given as X = Xez where X is angular velocity of rotation of impellers. The particle dynamic type simulation is performed to calculate the motion of particles in the liquid. The contact force as described in the preceding section is used to capture the major features of grain interactions. A particle dynamic model is employed to account for particle contacts of finite duration, in which the viscoelastic behavior of the particles is represented using a nonlinear Hertzian type model as presented in the preceding section. Using of a nonlinear viscoelastic model is essential for capturing granular structure formation and caking phenomenon in the mixing tank. Roughly 2 · 106 identical, slightly overlapping, spherical sodium formate particles with an initial diameter of dp0 = 1000 lm are used to fill the tank. The initial configuration of sodium formate particles in the tank is illustrated in Fig. 9(a). The free surface, which was a heap type before mixing begins, whose pick was located at h*  0.25 where h*, dimensionless height, is defined as h/H. h is height of liquid in the tank and H is the tank height. Eqs. (16) and (17) are integrated using 4th and 5th order embedded formulas from Dormand and Prince [34] with Dt = 5 · 108 s. The calculation of drag force acting on a particle requires knowledge of the local averaged values of the fluid velocity components at the position of the particle in the Lagrangian grid. Due to the numerical solution method detailed below these variables are only known at discrete nodes in the domain. Hence, a mass weighted averaging technique [35] has to be employed for calculating the averaged quantities in the Lagrangian grid. Fig. 9(a) shows the initial configuration of sodium formate grains in the tank. The particles are polydisperes with the average particle size of 530 lm. Using the diffusion layer model presented in the preceding section the particles shrink with different dissolution rates. Here, the Sherwood 0:33 number is estimated as Sh ¼ ½2ð1  /s Þ þ 0:86Re0:55  using LES technique whose sample results are prep Sc sented in Fig. 4. It is assumed that the shape of grains is always spherical. In actual situations, it is likely that particles become non-spherical for which no reliable expression exists for drag force such as that presented in (15). It is assumed that Dsol = 106 m2/s. These complexities are neglected at this stage. Fig. 9(b) represents an instantaneous configuration of the sodium formate particles after nearly t  1 s of simulation in which the rotation speed of the Rushton turbine is set to X = 400 rpm. As it can be seen from Fig. 9(b) that the granular system becomes polydisperse with the finer grains are located closer to the Rushton turbine. The collision time of fine particles, as illustrated in Fig. 10, having small impact velocity could be significantly large, which may be even of order of 103 s. This observation suggests that adding the adhesive force between sodium formate grains the collision becomes completely plastic with no restitution period after the approaching period during a collision. In this light, exceeding a critical size for grains phenomena, such as caking, might be prevented. However, the role of the adhesive force between the grains in the development of large solid structures in a brine-mixing tank merits further investigation. It is likely that non-spherical structures will form in actual situations. These complexities are neglected at this stage. Finally, results of this simulation have been compared with the previous work. Our results are in good agreement with Hartmann [36]. 3.1. Fluid–structure interaction Mixing processes in a mixing tank appears to be of an interdisciplinary nature comprising those physical phenomena which involve significant mutual interaction among hydrodynamics and structural forces. The problem involves significant nonlinearities and viscous effects but nothing that could be termed a systematic theory has been developed to date.


S.M. Mousavi et al. / Communications in Nonlinear Science and Numerical Simulation 13 (2008) 380–399

Fig. 11. (a) Grid for the rotor and blades. (b–d) Deformations of the elastic rotor caused by rotating liquid–particle system. The rotor and blades are color coded using the local displacement magnitude where the red color represents significant displacement as large as 3 mm. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

The advantage of the proposed approach is that issues such as fluid–structure interaction occurs in the mixing tanks can be investigated. Note that the flow of solution-particle mixture (as illustrated in Fig. 10) causes

S.M. Mousavi et al. / Communications in Nonlinear Science and Numerical Simulation 13 (2008) 380–399


the deformation of the rotor. This deformation of the rotor, in turn, could change the boundary condition of the fluid problem significantly. Using an approach outlined in [36] the deformation of the rotor and blades of the Rushton turbine in a mixing tank can be predicted. In this case, the analysis should provide a strong coupling between the dynamics of liquid–particle flow and the dynamics of structures. Using a finite volume approach for the Rushton turbine [36] a coupled system is developed in which the interaction forces between fluid and structure are accounted for and the resultant motions of the rotor is predicted as illustrated in Fig. 11. Presently, no details such as flow dynamics of sodium formate grains in a mixing tank at different rates of rotation of turbine can be provided using even the simplified approach described above. Limitations in the available computational resources urge using an alternative approach such as a continuum approach, which has been proposed in an earlier work [37]. Indeed, the discrete model developed in this paper provides substantial insight in the physics of a mixing tank, which could lead to improved constitutive equations for continuum models. 4. Conclusions In this work the development of a comprehensive particle dynamics model was presented for studying the dissolution of sodium formate grains in a mixing tank with Rushton turbine. A detailed overview has been provided of processes in liquid–particle flows to obtain insight into the hydrodynamics of a mixing tank and to arrive at a simplified but realistic model for a mixing tank. The simplified model includes an equation of motion for spherical sodium formate particles having repulsive and dissipative interactions with its neighbors. In addition, the liquid–particle interactions are taken into account using an empirical formula for drag force exerted by the liquid on the particles. In addition, the simplified diffusion layer model presented in order to estimate the dissolution rate and the concentration profile around a grain of sodium formate in a brine-mixing tank. The hydrodynamic model of the liquid phase is based on generalized Navier–Stokes equations for a liquid interacting with a solid phase. The governing equations of the gas phase also were solved using LES technique. Acknowledgement S.M. Mousavi wishes to acknowledge the support of the Finnish Academy (Award No. 211530).  2 Þ1=2 expðcjxi  ni j2 =D  2 Þ with filter width Appendix A. By applying the Gaussian filter Gðxi  ni Þ ¼ ðc=pD Di(i = x, y, z) (taken equal to the grid spacing) to Eqs. (24) and (25), the filtered equations are obtained given as ðqsol ð1  /s Þui Þ;i  qsol ð1  /s Þ;t ¼ 0;


ðqsol ð1  /s Þui Þ;t þ ðqsol ð1  /s Þui uj Þ;j ¼ ðð1  /s ÞpÞ;i þ ðð1  /s Þrij Þ;j  ð1  /s Þgqsol diz þ fi :


Any quantity such as ui in the flow domain can be decomposed as ui ¼ ui þ u0i , where u0i is the subgrid-scale part that accounts for the scales not resolved by the computational grid. Eqs. (24) and (25) look like the equation of motion of a fluid of variable density q = qsol(1  /s). Hence, as an alternative to the aforementioned approach discussed in Section 2, a Favre filter can be defined as ~ ui ¼

qui : q


This give rise to the alternative decomposition ui ¼ u~i þ u00i , where u00i represents the subgrid-scale part of ui based on Favre filtering. The Favre filtered governing equations are given as ;t ¼ 0; ðq~ ui Þ;i  q


ui Þ;t þ ðq~ ui ~ pÞ;i þ ðð1  /s Þ~ rij Þ;j þ sij;j  qgdiz þ fi : ðq~ uj Þ;j ¼ ðð1  /s Þ~



S.M. Mousavi et al. / Communications in Nonlinear Science and Numerical Simulation 13 (2008) 380–399

In (A.5), the molecular stress tensor, and the subgrid-scale stress tensor are given as 2 rij ¼  lsol up;p dij þ lsol ðui;j þ uj;i Þ 3 00 00 00 00 g ug ui ~ sij ¼ qð ~ uj  ~ uj þ ug uj þ ~ ug i~ i uj þ ui uj Þ: i~

ðA:6Þ ðA:7Þ

As a tentative first attempt at closure, the LES equations Eqs. (A.4) and (A.5) may be closed with a subgridscale model for, sij, given as [38]   1 1 1e 2 2 e e 2 e ~ ~ ~ ~ ~ S sij ¼ qð ~ ug  u Þ þ 2C qD ð S Þ  d S mn e ðA:8Þ u S S S mn Þ2 dij ; u  C 2 qD2i ð e i j i j 1 mn mn ij pp ij i 3 3 where C1 and C2 are dimensionless constants whose values may be determined by correlating the results of analysis such as that detailed in [39]. However, a highly resolved flow field is required to estimate the aforementioned constants accurately. Speziale et al. [38] found that the values of constants C1 and C2 for compressible turbulence are 0.012 and 0.0066, respectively. References [1] Mersmann A. Crystallization technology handbook. 2nd ed. New York: Marcel Dekker; 2001. [2] Armenante PM, Tsai D. Agitation requirements for complete dispersion of emulsions. Presented at the AICHE Annual Meeting, Washington, DC, November, 1988. [3] Bryant J. Characterization of mixing in fermenters. Adv Biochem Eng 1977;5:101–23. [4] Baldyga J, Bourne JR. A fluid mechanical approach to turbulent mixing and chemical reaction. Chem Eng Commun 1984;28:231–78. [5] Baldi G, Conti R, Gianetto A. Concentration profile for solids suspended in a continuous agitated reactor. AICHE J 1981;27:1017–20. [6] Barresi A, Baldi G. Solid dispersion in an agitated vessel. Chem Eng Sci 1987;42:2949–56. [7] Mactaggart RS, Nasr-El-Din HA, Masliyah JH. Sample withdrawal from a slurry mixing tank. Chem Eng Sci 1993;48:921–31. [8] Nasr-El-Din HA, Shook CA, Colwell J. A conductivity probe for measuring local concentration in slurry systems. Int J Multiphase Flow 1987;13:365–78. [9] Chudacek MW. Relationship between solid suspension criteria, mechanism of suspension, tank geometry, and scale-up parameters in stirred tank. Ind Eng Chem Fundam 1986;25:391–401. [10] Zwietering TN. Suspending of solid particles in liquid by agitators. Chem Eng Sci 1958;8:244–53. [11] Tatterson GB. Scaleup and design of industrial mixing processes. New York: McGraw-Hill; 1994. [12] Edward LP, Victor AAO, Kresta SM. Handbook of industrial mixing: science and practice. New Jersey: John Wiley & sons; 2004. [13] Jean M, Julian S. An introduction to turbulent flow. USA: Cambridge University Press; 2000. [14] Pope Stephan B. Turbulent flow. Cambridge, UK: Cambridge University Press; 2000. [15] Revstedt J, Fuchs L, Tragardh C. Large eddy simulation of the turbulent flow in a stirred tank. Chem Eng Sci 1998;53(24):4041–53. [16] Derksen J, Van den Akker HEA. Large eddy simulation on the flow driven by a Rushton turbine. AIChE J 1999;45(2):209–21. [17] Lu Z, Liao Y, Qian D, McLaughlin JB, Derksen JJ, Kontomaris K. Large eddy simulations of a stirred tank using the lattice Boltzmann method on a nonuniform grid. J Comput Phys 2002;181:675–704. [18] Moin P. Advances in large eddy simulation methodology for complex flows. Int J Heat Fluid Flow 2002;23:710–20. [19] Yeoh SL, Papadakis G, Lee KC, Yianneskis M. Large eddy simulation of turbulent flow in Rushton impeller stirred reactor with a sliding-deforming mesh methodology. Chem Eng Tech 2004;27(3):257–63. [20] Julsing HG, McCrindle RI. The recovery of precious metals from acidic effluents using sodium formate. Water Sci Technol 2000;42:63–9. [21] Groschuff E, Ber Dtsch. Chem Ges 1903;36:1783. 4351. [22] Sidgwick NV, Gentle JAHR. J Chem Soc 1922;121:1837. [23] Zamankhan P, Bordbar MH. Complex flow dynamics in dense granular flows—Part 1: Experimentation. J Appl Mech, in press. [24] Hirt CW, Nichols BD. Volume of fluid (VOF) method for the dynamics of free boundaries. J Comput Phys 1981;39:201. [25] Magnaudet J, Eames I. The dynamics of high Reynolds number bubbles in inhomogeneous flows. Annu Rev Fluid Mech 2000;32:659–708. [26] Menon S, Kim WW. High Reynolds number flow simulations using the localized dynamic subgrid-scale model. In: 34th aerospace science meeting, AIAA Paper 96-0425, Reno, 1996. [27] Chen L, Garimella SV, Reizes JA, Leonardi E. The development of a bubble rising in a viscous liquid. J Fluid Mech 1999;387:61–9. [28] Yamamoto Y, Potthoff M, Tanaka T, Kajishima T, Tsuji Y. Large-eddy simulation of turbulence gas–particle flow in a vertical channel: effect of considering inter-particle collisions. J Fluid Mech 2001;442:303–34. [29] Phan-Thien N. Understanding viscoelasticity. Heidelberg, Germany: Springer; 2002. [30] Silbert LE, Ertas D, Grest GS, Hasley TS, Levin D, Plimpton SJ. Phys Rev E 2001;64:051302. [31] Brilliantov NV, Spahn F, Hertzsch JM. Model for collisions in granular gases. Phys Rev E 1996;53(5):5382–92.

S.M. Mousavi et al. / Communications in Nonlinear Science and Numerical Simulation 13 (2008) 380–399


[32] Zamankhan P, Malinen P, Lepoma¨ki H. Application of neural networks to mass-transfer predictions in a fast fluidized bed of fine solids. AICHE J 1997;43(7):1684–90. [33] Cussler EL. Diffusion mass transfer in fluid system. New York, USA: Cambridge University Press; 1997. [34] Dormand JR, Prince PJ. A family of embedded Runge–Kutta formulae. J Comput Appl Math 1980;6:19–26. [35] Zamankhan P. Kinetic theory of multicomponent dense mixtures of slightly inelastic spherical particles. Phys Rev E 1995;52:4877. [36] Hartmann H. Detailed simulations of liquid and solid–liquid mixing. PhD thesis, Delft University, 2005. [37] Bordbar MH, Zamankhan P. Dynamical states of bubbling in vertically vibrated granular materials, Part II: Theoretical analysis and simulations. Commun Nonlinear Sci Numer Simul, in press, doi:10.1016/j.cnsns.2005.03.008. [38] Speziale CG, Erlebacher G, Zang TA, Hussaini MY. The subgrid-scale modeling of compressible turbulence. Phys Fluids 1988;31:940–2. [39] Bordbar MH, Zamankhan P. Dynamical states of bubbling in vertically vibrated granular materials, Part I: Collective processes. Commun Nonlinear Sci Numer Simul, in press, doi:10.1016/j.cnsns.2005.04.001.