- Email: [email protected]

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

a,b,*

, P. Zamankhan b, A. Jafari

b

a

b

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 ﬂuid in challenging environments such as the Barents Sea. In this paper, large-eddy simulations of a turbulent ﬂow in a solid–liquid baﬄed 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 ﬂow 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-ﬂowing 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

381

mixing tanks are required to fulﬁll 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 ﬂows 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 ﬂuid 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 signiﬁcant improvement on mixing tank design can be obtained by developing numerical models, which take into account the real ﬂow ﬁeld inside the vessel. Eﬀorts have been made during the last years towards the development of predictive methods based on computational ﬂuid dynamics (CFD) and capable of providing detailed information on ﬂow and turbulence ﬁelds of mixing tanks [12]. Suﬃciently 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 ﬂow scales in industrial applications, within a reasonable amount of time and memory space. A direct numerical simulation (DNS) of the turbulent ﬂow 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 ﬂow. 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 ﬂow dynamics in a mixing tank correctly. Recent applications of LES to the hydrodynamic simulation of industrial equipment such as baﬄed stirred tanks are presented in [15–19]. Large-eddy simulation (LES) is a simulation technique [14] based on the decomposition of the ﬂow ﬁeld 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 eﬀects 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 baﬄed 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) Diﬀerent views of the 6 impellers Rushton turbine with its characteristic dimensions.

382

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 ﬂuids. It is often used in challenging environments such as the Barents Sea. Formate brine solutions including those made from sodium are used in the oilﬁeld as drilling ﬂuids where ﬂuid viscosity must be low. These ﬂuids are made up of the formate solution itself and water soluble polymers. The ﬂuids 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 hydrosulﬁte, a common reductive bleaching chemical. More applications of sodium formate such as that in the recovery of precious metals from acidic eﬄuents 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 diﬃcult one because the individual particles of powder can stick together during storage and fuse into a solid mass that is diﬃcult to break up and hydrate in a brine-mixing tank. Published data for solubilities of sodium formate in water include Groscuﬀ [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 ﬂows useful for developing a simpliﬁed model for particle dynamic type simulations of dense liquid–particle ﬂows. 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 diﬀerence with the surrounding liquid. A simulation of gas–particle ﬂow is performed for predicting the ﬂow 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

383

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 ﬂow 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 ﬂow 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 ﬁnding an universal model is much higher when only small eddies are modeled. In the LES approach, the governing equations are obtained by spatially ﬁltering 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 eﬀort is independent of the Reynolds number if the small scales obey the inertial range spectrum and the near wall eﬀects are not important. 2.2. Governing equations Turbulence in a steadily rotating ﬂuid, such as found in brine-mixing tanks, plays an important role in solid dispersion. As illustrated in Fig. 3, suﬃciently 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 eﬀects of rotation and if it is suﬃciently 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 signiﬁcantly during a rotation period. 2.2.1. Air–sodium formate solution ﬂows 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-ﬂuid (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 ﬂows of air and solution as immiscible ﬂuids with interfaces eﬀects 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 eﬀects 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 ﬂow 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 ﬂows in a brine-mixing tank.

384

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

Fig. 3. (a) The instantaneous velocity vector ﬁeld 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 (deﬁned 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 ﬁeld in a cutting xy-plane lactated below of the impellers.

The ﬁltered 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 ﬁlter. 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 ﬁlter 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 ﬁlter. To remedy this further decomposition is required ui ui ui ui uj ¼ ð ui uj uj Þ þ uj ¼ Lij þ uj ;

ð3Þ

where Lij indicates interactions among the large scales. Using Eq. (3), the ﬁltered 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

385

In order to model the subgrid-scale tensor, Menon and Kim [26], suggested that based on the similarity assumption of diﬀerent scale sij = CkLij, a kinetic energy ktest can be deﬁned 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 diﬀusion 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 simpliﬁes 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 ;

ð9Þ

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

ð10Þ

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 ﬁnite 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 coeﬃcient and j represents the curvature of the interface and is deﬁned 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.

386

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 ﬂow 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 baﬄes. 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 baﬄes, 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 diﬀerencing scheme is used for spatial discretization and the time is advanced via a second-order implicit scheme. 2.2.2. Particle–solution ﬂows 2.2.2.1. 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

ð12Þ

where V pi is particle velocity, F fi is drag force from ﬂuid, 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 deﬁne 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 Þ

ð13Þ

2

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

N

ð14Þ

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 coeﬃcient, CLR is lift coeﬃcient 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

387

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

Properties

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 coeﬃcient Poisson’s ratio Instantaneous shear modulus Long time shear modulus The relaxation time Elastic modulus Chemical formula Molecular weight

Water

Temperature Kinematic viscosity Density

20 C 106 m2/s 998.2 kg/m3

Air

Temperature Kinematic viscosity Density

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

Steel

Poisson’s ratio Elastic modulus Density

0.29 2.1 · 1011 Pa 7890 kg/m3

The sodium formate particles have surface roughness on many diﬀerent 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 ﬁnite 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 coeﬃcient 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 ﬁrst 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 deﬁnes 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 deﬁned as es;jp s;jp s;jp Coulomb friction law describes friction between two colliding grains with a surface friction coeﬃcient, 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

388

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

Fig. 4. The complex velocity vector ﬁeld 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 diﬀerent 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 ﬁgure legend, the reader is referred to the web version of this article.)

PN torque on particle j may be deﬁned 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

ð17Þ

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 coeﬃcient, lp, may be expressed in terms of mesoscopic parameters. However, in the present attempt the friction coeﬃcient 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 coeﬃcient of friction is presented in the following. To obtain an expression appropriate to spherical particle-ﬂat 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

389

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 conﬁgurations in (a) and (b) are each separated by t 106 s.

2.2.2.2. Coeﬃcient of friction. This section addresses spills of sodium formats onto a ﬂat 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 ﬁnite duration. In this case, a model for collision between particles was used by artiﬁcially 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 eﬀect 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 diﬀerent 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.

390

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 ﬂat 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 conﬁguration 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 aﬀect the spill geometry. Fig. 6(b) represents the resulting shape of spill of glass beads using simulations where a static coeﬃcient of friction, lp = 0.65, is used. In order to obtain an angle of repose of hf 33, the coeﬃcient of static friction is set to 0.65 in the present simulations. A slightly higher coeﬃcient of static friction has to be used to model the sliding movements of particles over a ﬂat surface made of steel. Here, the coeﬃcient of static friction is set to 0.8 for grain steel sliding contacts. The base of the ﬂat surface is supported so that no movement is possible in any direction and the eﬀects 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). 2.2.2.3. Solubility. Solubility is the amount of sodium formate that will dissolve in the water. Solubility of sodium formate has special signiﬁcance 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 ﬁnd both the dissolution rate and the concentration proﬁle 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

391

Fig. 7. A diﬀusion model for the dissolution of a single grain of sodium formate in water. The concentration ﬁeld 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

ð18Þ

where r represents the distance from center of particle, and Cl is concentration of sodium formate within the boundary layer deﬁned 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 coeﬃcient of diﬀusivity of sodium formate (liquid) in water, Dsol, changes insigniﬁcantly 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 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ð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 ﬂux may be given as J ðtÞ ﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ Dsol =ptðC l0 C 1 Þ, whereas at long times, namely t R2p =pDsol , it simpliﬁes to a value, which depends on the grain size, as given by J ðRp Þ ﬃ Dsol ðC l0 C 1 Þ=Rp . This suggests that using the smaller size grains may increase the interfacial ﬂux. Here Cl0 is concentration of sodium formate at the solid surface, and C1 shows average concentration of sodium formate deﬁned 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 ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 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 pﬃﬃﬃ 2 2 2ðK sphere Þ ½1 pK sphere eðK sphere Þ erfcðK sphere Þ ¼ S; ð21Þ where S is dimensionless saturation deﬁned as facial speed given as pﬃﬃﬃﬃﬃﬃﬃﬃ dRp j ¼ ðK sphere Dsol Þt1=2 ¼ pﬃ : dt t

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

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

392

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

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

ð23Þ

where dBL is diﬀusion 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 diﬀusion layer thickness is an important factor in controlling the shape of particle dissolution proﬁles. 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 coeﬃcient may be estimated using Sh ¼ 2kRp =Dsol ¼ ½2ð1 /s Þ þ AReap Scb [32] where Sc is Schmidt number. Currently, no information exists for the coeﬃcient 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 simpliﬁed diﬀusion layer model presented in this section provides an approximation for estimating the dissolution rate and the concentration proﬁle around the grain of sodium formate in a brine-mixing tank. 2.2.2.4. Liquid phase. To predict the solution ﬂow ﬁeld 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;

ð24Þ

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

ð25Þ

The last term on the right hand side of Eq. (25) represents the particle eﬀect 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 diﬀusion layer thickness with saturation. Inset: The typical shape of particle dissolution proﬁle. 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

393

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

The governing equations of the liquid phase are solved using LES technique. The ﬁltered equations are similar to Eqs. (24) and (25) with the diﬀerence 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 ﬂow. In the present study, the model (A.8) is tested to clarify its relevance for dense water-sodium formate ﬂows in the mixing tank, for which the ﬂuctuating ﬁelds might not be small compared with those of the mean ﬁelds. 2.3. Numerical method A sliding mesh method has been used by which ﬂow 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 baﬄes. The unstructured grid in the impeller zone rotates with the impeller, whereas the grid in the tank remains stationary. The aforementioned grids slide past

394

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

Fig. 10. Diﬀerent views of a typical conﬁguration 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 baﬄes, 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 satisﬁed, 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 ﬁrst few iterations. It is worth noting that reﬁnement of the grids did not produce any signiﬁcant diﬀerences in the results. The ﬂow 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 speciﬁcation 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 ﬁnite 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 ﬁeld 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 diﬀerential equations (24), (25) and (A.8).

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

395

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 ﬂows 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 ﬂow 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 ﬁnite 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 ﬁll the tank. The initial conﬁguration 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 deﬁned 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 ﬂuid 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 conﬁguration of sodium formate grains in the tank. The particles are polydisperes with the average particle size of 530 lm. Using the diﬀusion layer model presented in the preceding section the particles shrink with diﬀerent 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 conﬁguration 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 ﬁner grains are located closer to the Rushton turbine. The collision time of ﬁne particles, as illustrated in Fig. 10, having small impact velocity could be signiﬁcantly 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 signiﬁcant mutual interaction among hydrodynamics and structural forces. The problem involves signiﬁcant nonlinearities and viscous eﬀects but nothing that could be termed a systematic theory has been developed to date.

396

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 signiﬁcant displacement as large as 3 mm. (For interpretation of the references to color in this ﬁgure legend, the reader is referred to the web version of this article.)

The advantage of the proposed approach is that issues such as ﬂuid–structure interaction occurs in the mixing tanks can be investigated. Note that the ﬂow 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

397

the deformation of the rotor. This deformation of the rotor, in turn, could change the boundary condition of the ﬂuid problem signiﬁcantly. 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 ﬂow and the dynamics of structures. Using a ﬁnite volume approach for the Rushton turbine [36] a coupled system is developed in which the interaction forces between ﬂuid and structure are accounted for and the resultant motions of the rotor is predicted as illustrated in Fig. 11. Presently, no details such as ﬂow dynamics of sodium formate grains in a mixing tank at diﬀerent rates of rotation of turbine can be provided using even the simpliﬁed 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 ﬂows to obtain insight into the hydrodynamics of a mixing tank and to arrive at a simpliﬁed but realistic model for a mixing tank. The simpliﬁed 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 simpliﬁed diﬀusion layer model presented in order to estimate the dissolution rate and the concentration proﬁle 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 ﬁlter width Appendix A. By applying the Gaussian ﬁlter Gðxi ni Þ ¼ ðc=pD Di(i = x, y, z) (taken equal to the grid spacing) to Eqs. (24) and (25), the ﬁltered equations are obtained given as ðqsol ð1 /s Þui Þ;i qsol ð1 /s Þ;t ¼ 0;

ðA:1Þ

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

ðA:2Þ

Any quantity such as ui in the ﬂow 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 ﬂuid of variable density q = qsol(1 /s). Hence, as an alternative to the aforementioned approach discussed in Section 2, a Favre ﬁlter can be deﬁned as ~ ui ¼

qui : q

ðA:3Þ

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

ðA:4Þ

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

ðA:5Þ

398

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 ﬁrst 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 ﬂow ﬁeld 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 ﬂuid mechanical approach to turbulent mixing and chemical reaction. Chem Eng Commun 1984;28:231–78. [5] Baldi G, Conti R, Gianetto A. Concentration proﬁle 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 ﬂow. USA: Cambridge University Press; 2000. [14] Pope Stephan B. Turbulent ﬂow. Cambridge, UK: Cambridge University Press; 2000. [15] Revstedt J, Fuchs L, Tragardh C. Large eddy simulation of the turbulent ﬂow in a stirred tank. Chem Eng Sci 1998;53(24):4041–53. [16] Derksen J, Van den Akker HEA. Large eddy simulation on the ﬂow 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 ﬂows. Int J Heat Fluid Flow 2002;23:710–20. [19] Yeoh SL, Papadakis G, Lee KC, Yianneskis M. Large eddy simulation of turbulent ﬂow 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 eﬄuents using sodium formate. Water Sci Technol 2000;42:63–9. [21] Groschuﬀ 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 ﬂow dynamics in dense granular ﬂows—Part 1: Experimentation. J Appl Mech, in press. [24] Hirt CW, Nichols BD. Volume of ﬂuid (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 ﬂows. Annu Rev Fluid Mech 2000;32:659–708. [26] Menon S, Kim WW. High Reynolds number ﬂow 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, Potthoﬀ M, Tanaka T, Kajishima T, Tsuji Y. Large-eddy simulation of turbulence gas–particle ﬂow in a vertical channel: eﬀect 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

399

[32] Zamankhan P, Malinen P, Lepoma¨ki H. Application of neural networks to mass-transfer predictions in a fast ﬂuidized bed of ﬁne solids. AICHE J 1997;43(7):1684–90. [33] Cussler EL. Diﬀusion mass transfer in ﬂuid 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.