Toward A Truly Multiscale Computational Strategy For Simulating ...

8 downloads 5250 Views 5MB Size Report
Oct 8, 2010 - Such a DNS is carried out in a periodic box, a dedicated forcing technique being used to impose the turbulent-flow conditions pertinent to a ...
10780

Ind. Eng. Chem. Res. 2010, 49, 10780–10797

Toward A Truly Multiscale Computational Strategy For Simulating Turbulent Two-Phase Flow Processes Harry E. A. Van den Akker* Multi-Scale Physics Department, Delft UniVersity of Technology, Delft, Netherlands

In chemical reactor engineering, simple concepts are used for describing the flow in a reactor. Turbulent two-phase flow processes, however, are characterized by a broad range of time and length scales. Two-phase reactors operated in the turbulent regime therefore qualify for multiscale modeling. Multiscale models are being developed in such diverging fields as chemical reaction engineering (among which are packed bed reactors, fluidized beds, and risers), chemical vapor deposition reactor modeling, turbulent single-phase and two-phase flow simulations (among which is combustion), and materials science. The common aspects of these multiscale models are highlighted: they all comprise a coarse-grained simulation for the macroscale and some type of fully resolved microscale simulation. Different simulation techniques are used however. Processes involving single-phase flow may require one of three computational fluid dynamics (CFD) techniques: direct numerical simulations (DNSs), large eddy simulations (LESs), and Reynolds averaged Navier-Stokes (RANS)-based simulations. For two-phase flows, two CFD options are open: Euler-Lagrangian (or particle tracking) and Euler-Euler (or two-fluid). The characteristics of all these approaches are discussed. One of the more interesting options in dealing with turbulent two-phase, i.e., multiscale, flow reactors is to run a DNS for the local small-scale processes. Such a DNS is carried out in a periodic box, a dedicated forcing technique being used to impose the turbulent-flow conditions pertinent to a specific position in the macro domain. Several such successful DNSs are reviewed, which all exploit lattice Boltzmann (LB) techniques. Some new and promising results from LB simulations for gas-liquid flow systems are presented. Finally, a truly multiscale simulation strategy is presented for turbulent two-phase flow reactors, which combines a coarse-grained simulation for the macro domain concurrently run with several DNSs for properly chosen positions in the domain. LB techniques are recommended. The crucial step of this strategysfeeding the results of the local DNS frequently back into the coarse-grained simulations until convergence is reachedsis described but has not been implemented yet. Introduction Chemical engineers usually design two-phase flow reactors on the basis of an idealized concept: steady, uniform, cocurrent or countercurrent, one-dimensional flow of the two phases, a homogeneous phase distribution across the cross-section of the reactor and along the height of the reactor, and a uniform and usually constant bubble, droplet or particle size. This is a gross simplification of the real flow and transport processes in twophase flow reactors. The traditional cures for dealing with nonideal flow characteristics comprise residence or age distributions, tanks-in-series models, compartment models, and axial and radial dispersion models (see, e.g., Levenspiel,1 Villermaux,2 Godbole and Shah,3 Mills and Chaudhari,4 and Van Geem et al.5). The same can be said about the two-phase model of a bubbling fluidized bed6 in which the two phases are the dense phase of the fluidized particles and the fast rising gas bubbles, respectively. These remedies, however, reflect the moderate capabilities of the mathematical tools available to chemical engineers in the preceding century. Building simple models reflecting the essentials of a process with a minimum of mathematics even grew into the trademark of the profession of chemical engineering. These days, prompted by both economic and ecological drivers, society urges for a more precise and more reliable description of the processes in the chemical process industry and the energy sector. As far as chemical reactors are concerned, the ideal has become to let the molecules meet at the right place at the right moment in time under the most favorable reaction * E-mail: [email protected].

conditions, to avoid a waste of energy and raw materials and the production of unwanted and useless or harmful byproducts, among which is CO2. In the meantime, computational power has grown tremendously, and computational fluid dynamics (CFD) has developed into a mature toolkit accepted and exploited by chemical engineers.7,8 This allows for a much more realistic treatment of the 3-D and transient character of the flow and transport phenomena in many chemical reactors and other process equipment.9,10 In addition, our understanding of the physics of multiphase flows has progressed substantially, not in the least thanks to the exploitation of advanced experimental techniques such as particle image velocimetry,11,12 laser-doppler anemometry,13,14 and optical bubble probes.15-18 Multiscale Character of Two-Phase Flows As a first result, chemical engineers have realized that, certainly under operating conditions usual in industry, onedimensional uniform multiphase flow does not exist at all. Uniform flow may not be stable against disturbances under most practical circumstances. Very special arrangements may be required to create and maintain steady, or smooth, uniform flow conditions in bubble columns.19 In reality, gas-liquid bubbly flows are typified by a wealth of transient coherent structures exhibiting varying degrees of vorticity.11,13,16,20,21 Van den Akker22 explicitly drew an analogy with single-phase buoyancy driven flow phenomena. These rather large (“mesoscale”) vortical structures or eddies show up as strong peaks at the low frequency end of a turbulent energy spectrum. At the high

10.1021/ie1006382  2010 American Chemical Society Published on Web 10/08/2010

Ind. Eng. Chem. Res., Vol. 49, No. 21, 2010

Figure 1. Turbulent energy spectra (i.e., the auto power spectral density function, APSD, in m2/s, as a function of frequency, in Hz) for singlephase flow (black), bubbly flow in a bubble column (blue), and bubbly flow in a transport line (red). Reprinted with permission from ref 23. Copyright 2001. Cambridge University Press.

frequency end of the spectrum, this spectrum sharply falls off, since at these short wavelengths, the turbulent kinetic energy is dissipated quite effectively into thermal energy. The turbulent energy spectrum of bubbly flows strongly resembles that of turbulent single-phase flows: see Figure 1, from the work of Mudde and Saito.23 The inevitable conclusion is that one really needs dynamic 3-D CFD simulations to represent the essential features of bubbly flows. Neither steady-state nor dynamic 2-D simulations with circular symmetry suffice.24 In bubbling fluidized beds,20,25 in the freeboards of bubbling fluidized beds,22,26,27 in turbulent fluidized beds22,28 and in gas-solids risers,29 the model of steady and uniform two-phase flow has been demonstrated to be too simplistic as well: coherent vortical structures and clusters or striations of particles occur abundantly and ubiquitously. In a long series of papers, the Li and Kwauk school in Beijing has reported about multiphase flow simulations comprising an energy minimization principle (see, e.g., Li and Kwauk20), which, however, is beyond the scope of the present paper. In itself, a turbulence spectrum as measured in gas-liquid bubbly flows23 is indicative of a wide range of time and length scales being characteristic of multiphase flows. The time and length scales associated with the vortical structures may relate to the characteristic dimensions, or integral scale, of the containing process equipment.21 In single-phase flow, the smallest, or “inner”, time and length scales have been named after Kolmogorov30 and depend on the fluid’s kinematic viscosity, ν, and on the (local) dissipation rate of turbulent kinetic energy, ε, according to ηK ) (ν3 /ε)1/4

(1)

The ratio of Kolmogorov scales to integral scales decreases with increasing Reynolds number.30 Increasing the Reynolds number at constant fluid properties can be done by increasing the velocity scale and/or by increasing the length scale of the equipment. Increasing the velocity scale at constant equipment

10781

size requires increasing the power input (which in a steady state is equal to ε) and implies a higher turbulence level; the result is that the fine scales shift to lower values. Upscaling process equipment (from lab scale via pilot scale to full industrial scale) at constant specific power input (or constant ε) increases the Reynolds number and the ratio of Kolmogorov scales to integral scales as well. The wide range of scales from equipment size (or integral scales) down to the Kolmogorov scales is therefore an inherent characteristic of (single-phase) turbulent flows. In two-phase flow, the abundant presence of vortical coherent structures discussed above is an additional source of turbulence22,29 across a wide range of scales. Therefore, the whole notion of twophase flow turbulence clashes with the aforementioned classical models of uniform 1-D flow. As such, multiphase flow qualifies for multiscale analysis in terms of both time and length scales indeed. There is much debate in the literature as to how turbulence is created in two-phase flow. The role of the mesoscale coherent structures has already been highlighted above. Agrawal et al.29 investigated the role of mesoscale structures in rapid gas-solid flows, identified several mechanisms for cluster formation, and suggested that local disturbances (small clusters and streamers) are at the origin of turbulence and grow toward larger scales. In the context of gas-liquid bubbly flows, the discussion focuses on whether shear is the most important mechanism creating turbulence or whether the slip velocity of the bubbles, drops, or particles (or the Reynolds number based on this slip velocity) is the dominant factor in the production of turbulent stresses (see also Van den Akker31). Several authors report that the motion of bubbles damps the turbulence levels and may even promote a return to isotropy (see also Van den Akker31). In general, small (single) particles are found to attenuate turbulence while large particles augment turbulence,32 both effects increasing with increasing mass loading. An important parameter governing the particle-turbulence interaction seems to be the ratio of particle size to integral length scale of the turbulent carrier phase,32 although other authors, e.g., Sato et al.,33 suggest the ratio of particle size to Kolmogorov length scale to be relevant. In the literature, several reports can be found reporting on numerical simulations investigating the particle-turbulence interaction. A detailed discussion of these papers is beyond the scope of the present paper however. It suffices to refer to Ten Cate et al.34 who were capable of simulating the dual effect of solid particles on the turbulence spectrum by means of a direct numerical simulation in a small periodic box. As a matter of fact, the interaction between the spectrum of turbulent eddies on the one hand and the bubbles, drops or particles on the other hand is another reason for a multiscale analysis of multiphase flows. The abundant presence of vortical structures along with the concept of a turbulent energy spectrum raises the suggestion that these flows exhibit strongly random features. Rather than random, these flows may be classified as chaotic, however. Multiphase flow is governed fullysi.e. over the full range of time and length scalessby the familiar transport equations, or conservation laws, for mass, momentum, energy, angular momentum, et cetera. That is why CFD has turned into such a powerful tool for simulating multiphase flow phenomena and processes. In chemical reactors, the chemical reactions take place on the molecular scale where the chemical species have to meet one another. Molecular diffusion has to take care of this. In

10782

Ind. Eng. Chem. Res., Vol. 49, No. 21, 2010

gases, the time scale for diffusion L2/Dswhere D is the molecular diffusion coefficientsis more or less equal to the time scale for molecular momentum transport (due to viscosity) L2/ νswhere ν is the kinematic viscosity coefficient. In liquids, however, mass diffusion is much slower than momentum diffusion; as a result, diffusion of chemical species is governed by the Batchelor length scale ηB that is related to the Kolmogorov length scale according to ηB ) ηKSc-1/2

(2)

where Sc denotes the nondimensional Schmidt number, which is the ratio of kinematic viscosity to mass diffusion coefficient. For liquids, the Schmidt number is on the order of 102 to 103, and hence, the Batchelor scale is smaller than the Kolmogorov scale. In any liquid-phase system where molecular diffusion plays a role, among which are chemical reactors, the spectrum of relevant length scales may extend 1 or 2 orders of magnitude below the Kolmogorov scale. This fact increases the task of multiscale modeling in turbulent two-phase flows even further. Further Multiscale Processes Turbulent two-phase flows are not the only example of multiscale systems. In principle, any chemical reactor in which transport phenomena (mixing, diffusion of species, heat supply, heat withdrawal) affect yield and/or selectivity is a multiscale system. Usually, chemical engineers use dimensionless numbers such as the (second) Damko¨hler number to simplify the analysis of such reactors. Dimensionless numbers may be conceived as the ratio of time scales or of length scales, the second Damko¨hler number DaII denoting the ratio of the transfer time scale to the time scale of the chemical reaction. Deen’s textbook Analysis of Transport Phenomena35 provides an excellent introduction to several classical scaling and approximation techniques, such as the similarity method, regular and singular perturbation analyses, and the integral approximation method. Such techniques sharpen the insight into which mechanism may be dominant in complicated processes involving multiple effects and multiple time, length, and velocity scales and make it possible to obtain reasonable solutions for difficult problems. Such approximate techniques therefore have great educational value. These days, however, computational resources and algorithms have become widely availablesthereby eliminating the immediate need of these approximate analytical techniques. This development in its turn has triggered an increasing demand for high-fidelity predictive computer codes for design and engineering purposes. The leap in computational power and numerical algorithms since the 1990s allowed for the development of a multiscale approach in which the different scales in a process are modeled. Figure 2 (from the work of Lerou and Ng36) illustrates the wide ranges of time and length scales in chemical reaction engineering. Figure 3 (from the work of Li et al.37) illustrates the multiscale character of many processes of interest in a different way. Lerou and Ng36 make a plea for a multiscale approach for commercializing a process by combining results from models and experiments as to the processes taking place on the various scales from the plant scale down to the molecular scale. This synthesizing step was realized by incorporating knowledge and tools traditionally belonging to different fields of expertise such as catalysis and fluid mechanics.38 Prasad and Vlachos39 named this trend a paradigm shift toward using first-principles-based models, which in many cases are multiscale in nature. Vlachos

Figure 2. The ranges of length and time scales typical of a chemical process or plant. Reprinted with permission from ref 36. Copyright 1996. Elsevier.

et al.40 gave a report in terms of a multiscale simulation ladder, as a means to illustrate two-way information traffic between the various time and length scales, and explicitly mentioned CFD as part of the macroscopic scale modeling. A nice prototypical example of exploiting sophisticated computational (numerical) techniques based on first principles (i.e., the 3-D Navier-Stokes equation) for dealing with the complex interaction of turbulent fluid flow, transport phenomena (mixing, heat transport), and multiple chemical reactions may be found in a few papers on the potential occurrence of hot spots in tubular single-phase low-density polyethylene reactors.41,42 3-D turbulent flow and heat and mass transport phenomena in such a reactor were explicitly resolved for a substantial part of the eddy spectrum by means of a large eddy simulation. Models were still required, however, to take the subgrid scale processes into account: both the small-scale turbulent eddies and the micromixing of the chemical species within the grid cells need delicate and dedicated models. Their work illustrates how fluid mechanics has become an integral part of reaction engineering indeed, as foreseen by Wei38 as early as 1991 (cited by Lerou and Ng36). A relevant question is whether or not the computational simulation of two-phase flow systems could derive advantage from the progress made in multiscale modeling in other fields of chemical engineering. As a matter of fact, multiscale methods are explicitly being used these days in describing several more processes, such as: • heterogeneous catalytic systems36,40 • high-pressure polymerization reactors41,42 • packed bed reactor modeling43-45 • fluidized beds, circulating fluidized beds, and gas-solid risers29,46-50 • turbulent combustion processes51,52 • low-pressure processes, such as in chemical vapor deposition (CVD) reactors and micropropulsion systems53-56 • materials modeling (crystal growth, defect dynamics, surface diffusion, adsorption and desorption, and phase transitions)57 All of these examples are characterized by the occurrence of multiple relevant time and length scales. Figures 4 and 5 illustrate this for combustion processes and CVD reactors, respectively. As a matter of fact, the time scales of chemical reactions are much shorter than the time scales associated with flow, mixing, diffusion, and/or conduction. This implies that apart from processes taking place in the chemical regime,

Ind. Eng. Chem. Res., Vol. 49, No. 21, 2010

10783

Figure 3. A hierarchy of time and length scales in terms of structures, physical scales, and levels. Reprinted with permission from ref 37. Copyright 2005. The Institution of Chemical Engineers in Rugby, U.K.

Figure 4. The multiple scales playing a role in combustion. Moving clockwise: from furnace down to the molecular scale. Reprinted with permission from D. J. E. M. Roekaerts, Delft University of Technology.

multiple time and length scales play a role in most industrialscale chemical reactors. Also in physical phenomena in materials at the molecular and atomic scales, multiple time scales can be discerned:57 atoms or molecules vibrate around specific locations at which free energy is minimal and which are separated by high free energy barriers; the jumps crossing these barriers take place at long time intervals compared with the time scale of the thermal vibrations.

All of these prototypical processes are presumed to exhibit a large separation of time and/or length scales. That is why, in the classical chemical engineering approach, simplifications based on a proper analysis of time scales yields sufficiently reliable results.35 In the world of numerical analysis, such separation in scales is known as stiffness.52,54,57 The short time scale imposes time steps in the numerical treatment to be excessively small with respect to the rare events or slow

10784

Ind. Eng. Chem. Res., Vol. 49, No. 21, 2010

Figure 5. The multiple scales playing a role in thin film processing by chemical vapor deposition for chip manufacturing. Moving clockwise: from wafersize down to the molecular size. Reprinted with permission from C. R. Kleijn, Delft University of Technology.

transport processes. In computationally simulating the various processes listed above, different strategies are exploited to overcome the aforementioned stiffness. Something similar may be said as to the occurrence of multiple length scales: coarse graining is the term29,40,57 denoting a common remedy that typically results in reduced computational cost at the expense of a smaller degree of physical fidelity. In his review paper on multiscale methods in turbulent combustion, Echekki52 successively treats the various strategies used in this field for overcoming the limitations of resolving all the time and length scales. First of all, he discusses model adaptivity: techniques for modifying the governing equations to the effect that the simulation becomes feasible in a “reasonable” amount of time and fits the computational power available. These techniques decouple and eliminate (or erase) ranges of scales and are common practice, even outside the turbulent combustion field (see, e.g., Agrawal et al.29). This practice is the equivalent of coarse graining mentioned above. Typical results for single-phase flows are Reynolds averaged Navier-Stokes (RANS)-based simulations and large eddy simulations (LESs), which will be discussed in greater detail further on in this paper. A second class of multiscale strategies Echekki distinguishes comprises mesh adaptivity methods in which the governing equations are applied across all relevant scales of the flow. Adaptive multilevel (or multigrid) implementations based on a hierarchy of grid resolutions (also denoted as adaptive mesh refinement techniques) as well as spectral or wavelet-based (i.e., multiresolution) methods capable of switching between information in physical space and wavenumber space belong to this class. More details can be found in Echekki’s review paper.52 A third class according to Echekki exploits flame embedding methods which are not only very popular in the field of turbulent combustion but also very specific to just this field. They therefore are beyond the scope of this paper. Hybrid techniques, however, in which low dimensional models are embedded in the framework of a LES, represent a multiscale strategy worth highlighting. In these hybrid techniques, the processes taking place within a LES grid cell are simulated by a low-dimensional model. By doing so, a useful separation of scales is obtained,

which results in a substantial savings in computational time. Further on in the paper, we will come back to these methods. Another issue related to length scales occurs in the low pressure processes referred to above (CVD reactors and micropropulsion systems). There, the Knudsen number Knsdenoting the ratio of the mean free path of the individual molecules to the equipment dimensionsmay become order 1 or larger. As a result, a continuum description may break down locally,55,58 and one may have to switch to direct simulation Monte Carlo (DSMC) methods in order to reproduce the characteristics of these dilute systems. A typical example of a hybrid approach combining a Navier-Stokes (NS) solver with a DSMC method is the paper on the multiscale simulation of rarefied gas flows due to Abbate et al.56 In CVD reactors, operated under laminar-flow conditions, the usual goal is to deposit multiple, very thin, and uniform layers on wafers from precursors in the gas phase. In some of these process steps, submicrometer trenches and contact holes have to be filledspreferably completelyswith a particular material; the submicrometer size of these features, combined with the low pressures applied, may again require Monte Carlo (MC) methods.53,59 An important issue is the coupling between microscale (feature) and macroscale (reactor) processes (see also Figure 5): the boundary conditions for the microscale models are dictated by the local macroscale conditions (species concentration, temperature), while the microscale processes impose the boundary conditions for the macroscale simulation of flow and transport processes.55 It strongly depends, however, on the kind of phenomena and processes at stake which methods or equations are used for describing them. In materials modeling, various types and levels of MC algorithms have been in vogue since the 1970s. Chatterjee and Vlachos57 mention the kinetic MC method for simulating the events on a microscopic lattice, a probabilityweighted MC method for various diffusion processes, and a spatial coarse-grained MC method for simulations at larger length and time scales thanks to grouping lattice sites into coarse cells. As a result of the tempestuous growth in computational power since the 1990s, however, molecular dynamicssvery

Ind. Eng. Chem. Res., Vol. 49, No. 21, 2010

truly a first-principles toolshas gradually developed into the computational method of choice in the world of microscopic material modeling.57 On the contrary, in chemical reaction engineeringswith catalytic reactors, fluidized beds, and turbulent flames as typical instances of interestsmany fine-scale processes are embedded in a flowing fluid. To deal with this feature of flow processes, Kuipers and his group at Twente University46-48 developed a hierarchy of models for fluidized beds: these will be discussed further on in this paper in more detail. Their hierarchy of models is more or less analogous to the set of MC methods in materials modeling referred to above. Generally speaking, the motion of fluids is perfectly described by means of the classical Navier-Stokes (NS) equation. Although this equationsalso based on first principlesswas first written down in the 19th century, our understanding of it remains minimal, particularly because this soft-field vector equation is 3-D and strongly nonlinear. Dedicated software is needed to find accurate numerical solutions. Special versions of the NS equation apply to turbulent flows and multiphase flows. The unique character of the NS equation implies that just a limited advantage may be gained from the developments in multiscale modeling in other fields of chemical engineering (except combustion) which largely exploit MC techniques. The remainder of this paper will focus on the various options for dealing with the multiple time and length scales in chemical processes and physical operations carried out under turbulent two-phase flow conditions. Such flow conditions have a big impact on the circumstances (such as species concentrations and temperature) under which many chemical reactions or physical transformations proceed. The important issue is therefore how we can take the effects of turbulence into account when we like to mimic or predict the yield of chemical reactions and physical operations such as mixing and mass and heat transfer. It would be nice if the full range of scales present in a turbulent flow could be simulated in all detailsboth in singlephase flow and in two-phase flow reactors. CFD techniques, however, are still not capable of fully resolving a highly turbulent single-phase flow in most engineering equipment within a reasonable time span, let alone a turbulent two-phase flow. For the sake of a good comprehension of the options available, a short introduction to the various types of singlephase flow simulations is essential. The reader is also referred to Kuipers and Van Swaaij7 and to Van den Akker.8,10 Various Types of Single-Phase Flow Simulations For computationally simulating single-phase flows, three options are at our disposal: direct numerical simulations, which resolve all scales in the flow, Reynolds averaged NavierStokes-based simulations, which reproduce the mean flow characteristics only, and large eddy simulations, which result in a dynamic representation of a turbulent flow including a substantial part of the length and time scales in that flow. These three options will be summarized first. Direct Numerical Simulations (DNSs). In a direct numerical simulation (DNS), all scales of the velocity field v are simulated by means of the classical NS equation: ∂v 1 + v · ∇v ) - ∇p + ν∇2v ∂t F

(3)

in which p represents pressure and F and ν denote the density and viscosity of the fluid, respectively. Such a DNS is perfectly suited for two purposes:

10785

• on the macroscale, i.e., when the flow domain of interest comprises the whole piece of process equipment: e.g., a single solid spherical particle settling under the effect of gravity in a box filled with liquid originally at rest60,61 or the laminar flow field in a static mixer, even if the conditions are such that the flow is unsteady62 • on the microscale: when the interest is in resolving locally all details of the flow, e.g., around and between colliding particles immersed in a turbulent flow of a carrier phase (see Figure 7, from the work of Ten Cate et al.34); the flow domain is then kept restricted to a small box with periodic walls; this approach will be discussed further on in this paper in much greater detail For simulating strongly turbulent flows in a flow or process device, however, a full DNS resolving all eddies of the turbulent spectrum including the smallest would require the spacing of the computational grid to be on the order of the Kolmogorov length scale. Since the latter is as small as Re-3/4 times the macroscopic dimension of interest,30 a 3-D DNS of a turbulent flow would scale with Re9/4 (e.g., Wilcox63 and Moin and Kim64). This really limits the applicability of DNS to rather low Reynolds numbers. Sommerfeld and Decker65 report a DNS of the (marginally) turbulent flow in a lab-scale stirred vessel at a Reynolds number Re of 8000 to take approximately 3 months on eight processors with more than 17 GB of memory. Very recently, Gillissen and Van den Akker66 carried out a DNS of a stirred vessel at Re ) 7300 on a grid of some 2.9 × 109 grid cells; this took 128 CPUs and 336 GB of computer memory at the Dutch National Supercomputer SARA in Amsterdam. Such computational power is far beyond what chemical engineers in the industry can afford, while many industrial processes are carried out at (much) higher Reynolds numbers. Usually, a DNS is impractical for that reason. Reynolds Averaged Navier-Stokes (RANS) Simulations. Hence, turbulent flows are usually simulated with the help of the Reynolds aVeraged NaVier-Stokes (RANS) equation,30 which deliver an aVeraged representation of the flow only. Due to the averaging procedure, the NS equation comprises an additional term representing the extra contribution of the full spectrum of turbulent scales (eddies) to the momentum redistribution across the flow domain. Frisch67 recollects that as early as 1870 Boussinesq stressed that turbulence greatly increases viscosity (i.e., momentum exchange) and proposed the concept of a so-called eddy or turbulent viscosity, νt. The eventual equation for the averaged velocity field V including the averaged pressure field P looks like ∂V 1 + V · ∇V ) - ∇P + (ν + νt)∇2V ∂t F

(4)

In its turn, the turbulent viscosity may be position-dependent and requires modeling, for which a k-ε model is very usual: νt ) C µ

k2 ε

(5)

where k is the concentration of turbulent kinetic energy in J/kg (or m2/s2) and ε is the rate of dissipation (in W/kg, or m2/s2/s) of this turbulent kinetic energy, with Cµ denoting an empirical dimensionless coefficient which is to be treated as a universal constant. The two concentrations k and ε are generally interpreted as the most important parameters characterizing a turbulent flow field. In their turn, their spatial distributions within the turbulent

10786

Ind. Eng. Chem. Res., Vol. 49, No. 21, 2010

flow domain may be calculated from suitable transport equations.10,68 Various other closure models for taking the effect of the full spectrum of turbulent eddies on the average flow field into account are widely used, among which are so-called Reynolds stress models. Whatever choice is made for such modeling, it may lead to poor results, particularly with respect to small-scale phenomena, since many of the latter are nonlinearly dependent on the flow field.69 RANS-based simulations exploiting some k-ε models just reproduce the average flow field along with the spatial distributions of k and ε. Compared to the fine mesh used in a DNS and the resulting detailed field information, a RANS-based simulation may be named “coarse-grained” in the field of computational science. As such, RANS-based simulations are excellently suited for identifying dead zones, recirculatory flows, shortcircuiting between the entrance and exit, and further undesired flow features. In transient RANS-based simulations, however, it is not a priori clear which part of the fluctuations is temporally resolved and which part is taken care of by the turbulence model. This inherent property of RANS-based simulations especially raises concerns in the case of flows exhibiting no clear spectral separation between low-frequency coherent motions (such as macro-instabilities, precessing vortices, and trailing vortices) and turbulent fluctuations (making part of the cascade of eddies or whirls typical of turbulence). Large Eddy Simulations (LESs). The exponential increase in computational resources (along with a drastic drop in price for computer hardware) has led to the ability to solve industrial flows by means of large eddy simulations (LESs). The basic idea of LES is that the turbulent scales larger than the grid spacing ∆ are resolved explicitly, while eddies smaller than ∆ are filtered out and not resolved. This means that grid spacing acts as a low pass filter. These small eddies within each grid cell do contribute, however, to the redistribution of momentum into the adjacent grid cells and from there across the whole flow domain. The resulting NS equation to be solved for a LES now runs as ∂v˜ 1 + v˜ · ∇v˜ ) - ∇p˜ + ν∇2v˜ - ∇ · τ ∂t F

(6)

where v˜ and p˜ now denote the variables to be resolved at the detailed level of the grid spacing. The latter term in eq 6 represents the effect of the so-called subgrid scale (SGS) eddies. It is then assumed that due to this separation in scales, so-called subgrid scale (SGS) modeling is largely geometry-independent due to the presumed universal behavior of turbulence on the small scales. This subgrid scale contribution is modeled by means of an SGS eddy viscosity coefficient, νe. Usually, the model introduced by Smagorinsky70 suffices; it presumes that νe is proportional to ∆2 and to the magnitude S of the rate of strain in the grid cell of interest (see, e.g., Derksen71 and Van den Akker10). To take full advantage of the concept of LES, which is particularly aimed at resolving a great deal of the time and length scales of a turbulent flow field, fine grids should be used. This implies that LESs applied to a flow domain of some size are computationally quite demanding. Restricting LES to 2-D in view of saving computational time does not make sense physically, as the dynamics of turbulent flows in process equipment is inherently 3-D. In an attempt to keep simulation times restricted, it is not an option to take refuge in coarse grids either, as LES on a coarse grid does not make sense physically: the grid spacing should be so small that at least (a substantial) part of the inertial subrange can be resolved.

Compared to RANS-based simulations, LESs yield much more detailed simulation results indeed, not only because the grids used easily comprise millions of grid cellssthese numbers being larger than what is common in RANS simulations by at least 1 order of magnitudesbut also because the simulations are inherently transient and reproduce the dynamics of a large proportion of the wide spectrum of eddies. In terms of computational resources needed, LESs are positioned somewhere in between DNS and RANS-based simulations. Derksen72 carried out LES for a stirred vessel with Re ) 105 on a grid of 1.4 × 107 nodes, while a DNS would have required a grid of as many as 1012 nodes, the latter number being far beyond current computational capabilities. Comparisons with quantitative experimental data have increased the confidence in LES (see, e.g., Hartmann et al.).73 This does not apply to the hydrodynamics only but relates to the physical operations and chemical processes carried out in turbulent stirred vessels as well.10 After this review of the various CFD options for single-phase flows, it is time to present two essentially different ways of computationally simulating (turbulent) two-phase flows. Types of Two-Phase Flow Simulations Euler-Lagrangian Simulations. Often, dilute two-phase flows in which one phase is finely dispersed in a continuous carrier phase are computationally simulated by means of a Euler-Lagrange (E-L) method: first, the flow of the continuous phase is calculated in a Eulerian framework according to one of the three aforementioned options (RANS, LES, or DNS); then, the flow behavior of the second, dispersed phase is mimicked by tracking the path of individual point particles in a Lagrangian way. This can be done one-way (the particles feel the flow of the carrier phase, while the latter does not feel the particles74) or two-way and iteratively.75 The consequence of treating particles as point particles is that the detailed flow between the particles in response to the presence and motion of the particles remains unresolved. Without much discussion, one may anticipate that, in any case, particle inertia, gravity, and some fluid-particle interaction force need to be part of the equation of motion describing the motions and paths of the particles. Since the major fluid-particle force may be the drag force, the fluid’s velocity field is of primary importance, the turbulent velocity fluctuations inclusively, whether or not they are resolved in the simulation; fluctuations in pressure and stresses may be secondary. Supplying stochastic variations on top of the resolved velocity field mimics the unresolved fluctuations and brings the expected seemingly erratic paths of the particles about. Of course, the role of the artificially introduced stochastics for mimicking the effect of all eddies in a RANS-based particle tracking is different from and much more pronounced than that for mimicking the effect of just the SGS eddies in a LES-based tracking procedure. For more details on E-L simulations, the reader is referred to Patankar and Joseph76 and Van den Akker.10 Sommerfeld has exploited E-L widely: see, e.g., Sommerfeld and Decker.65 Very often, it is difficult to estimate a priori if more exotic fluid-particle forces such as the Saffman lift force, the Magnus force, and the history force may play a significant role either globally or locally. For a concise summary about these forces and for expressions for these forces, the reader is referred to, e.g., Derksen.72 Added mass may be important for bubbly flows. It is obvious that in such a Lagrangian approach distributions in particle size may easily be taken into accountsthough at the cost of computer time.

Ind. Eng. Chem. Res., Vol. 49, No. 21, 2010

For the hydrodynamics forces acting on the particles, mostly single-particle expressions are used; this implies that hydrodynamic interactions between particles are ignored completely. In many cases, direct interactions of particles owing to collisions are ignored as well. All thissalong with the computational burden that increases linearly with the number of particles being trackedsmay limit the practical applicability of the method to dilute systems with relatively low volume fractions of dispersed phase and/or to flow domains of small size. At the same time, the moving, usually solid, particles may be subjected to physical77 or chemical78 processes. For more dense flows, tracking all individual particles becomes computationally very intensive. The results reported by Chen et al.,79 who tracked 0.1 million particles in a gas flow, are very promising however. A cheaper way out may be the clustering of particles to groups, e.g., to mimic mass loading effects in gas-solid cyclones.75 The best alternative might be to pursue Euler-Euler (E-E), or two-fluid, simulations. Euler-Euler Simulations. In the complete Eulerian (E-E) description of multiphase flows, the dispersed phase may well be conceived as a second continuous phase that interpenetrates the real continuous phase, the carrier phase; this approach is often referred to as two-fluid formulation. The resulting “simultaneous” presence of two continua is taken into account by their respective volume fractions. All other variables such as velocities need to be averaged, in some way, in proportion to their presence; various techniques have been proposed to that purpose. In the two-fluid formulation, the motion or velocity field of each of the two continuous phases is described by its own continuity equation and NS equation (see e.g., Rietema and Van den Akker80 and Agrawal et al.29). The two NS equations look like ∂Vc + RFcVc · ∇Vc ) ∂t -R∇P + R(µ + µt)∇2Vc + RFcg - Fi RFc

(7)

∂Vd + (1 - R)FdVd · ∇Vd ) ∂t -(1 - R)∇P + (1 - R)(µ + µt)∇2Vc + (1 - R)Fdg + Fi (8) (1 - R)Fd

In both momentum balancesseqs 7 and 8sa phase interaction force Fi between the two interpenetrating continuous phases occurs predominantly, of course with opposite sign. Two-fluid models therefore belong to the class of two-way coupling approaches. The continuum formulation of the phase interaction force should reflect the same effects as experienced by the individual particles and discussed above in the context of the Lagrangian description of dispersed two-phase flow. One therefore has to decide here too which components of the phase interaction force (drag, virtual mass, Saffman lift, Magnus, history, stress gradients) are relevant and should be incorporated in the two sets of NS equations. Several authors, see, e.g., Oey et al.,81 have reported about the effects of ignoring certain components of the interaction force in the two-fluid approach. The question how to model in the two-fluid formulation (lateral) dispersion of bubbles, drops, and particles in swarms is relevant as well: various models are around. Another important issue in two-fluid models is about modeling the turbulent stresses under two-phase conditions.31,82 At this moment, there is still no consensus on a universal twophase turbulence model. Generally, turbulence in the continuous

10787

phase may be generated by shear due to large-scale velocity gradients felt by the continuous phase itself (just like in singlephase flows) as well as by the presence and relative motion of the dispersed phase particles. The ratio at which these two mechanisms contribute to the generation of turbulence may be an important factor in drafting a universal model. In addition, the dispersed phase may exhibit a turbulent flow behavior in response to the turbulent motions of the continuous phase in which it is embedded; this response may depend on several time scales and their interaction.81 Lance et al.83 suggested that the motion of bubbles promotes a return to isotropy (see also Van den Akker31). Drastic assumptions as to how to model the turbulent character of the two-phase flow are not uncommon. A universal model is not available right now. In addition, in two-fluid models, particle size is just a parameter in the phase interaction force for which an empirical expression is usedsusually restricted to the steady-state drag force for a single particle in an unbounded uniform 1-D flow. Which (average) particle size is appropriate and valid for the whole vessel or reactor of interest is up to the judgment of the person carrying out the simulation. This type of simulation does not allow for spatial variations in bubble or drop size in response to spatial variations in shear rate or dissipation rate of turbulent kinetic energy ε. In dense systems such as those encountered in solid suspensions, particle-particle interaction may be important as wellsthough ignored in eq 8. Then, the closure of solid-phase stresses is an important issue for which kinetic theory models and solid phase viscosity may be instrumental.84 In addition, fluidized beds may be dominated by particle-particle collisions (see further on). As a matter of fact, in comparison with the E-L approach, the complete Eulerian (or E-E) approach may better comply with denser two-phase flows, i.e., with higher volume fractions of the dispersed phase, when tracking individual particles is no longer doable in view of the computational times involved and the computer memory required, and when the physical interactions become too dominating to be ignored. Under these circumstances, the motion of individual particles may be overlooked, and it is wiser to opt for a more superficial, “coarsegrained” strategy that, however, still has to take the proper physics into account. A proper assessment of E-L versus E-E simulations was presented by Portela and Oliemans.85 Most E-E simulations are of the RANS type, although recently results obtained via an E-E, or two-fluid, LES have also been reported.82,86 Agrawal et al.29 concluded that coarsegrid simulations of gas-particle flows must include subgrid models, to account for the effects of the unresolved mesoscale structures. A similar conclusion was drawn by Wang et al.50 for circulating fluidized beds. In his review paper, Wang49 distinguished between six categories of subgrid scale models for fluidized beds. In the case of droplets and bubbles, particle size and particle number density may respond to variations in shear or energy dissipation rate. Such variations are abundantly present in turbulent stirred vessels. In fact, the explicit role of the revolving impeller is to produce small bubbles or drops, while in substantial parts of the vessel coalescence may be able to increase bubble or drop size again due to locally lower turbulence levels. Particle size distributions and their spatial variations are therefore commonplace and unavoidable in industrial mixing equipment. This seriously limits the applicability of common E-E models exploiting just a single value for particle size.

10788

Ind. Eng. Chem. Res., Vol. 49, No. 21, 2010

Some papers (e.g., Venneker et al.87) present results obtained by combining E-E with population balance models taking into account varying particle size distributions, but this combination is still quite cumbersome: it requires compromises about the flow field simulation and struggles with convergence. Most authors take refuge in some method of moments in view of allowing for a particle size distribution (e.g., Petitti et al.88), but even then the simulation is tedious and time-consuming. Another way out is to adopt a multifluid or multiphase approach in which various particle size classes are distinguished, each with its own NS equation, with mutual transition paths due to particle break-up and coalescence. A Hierarchy of Multiscale Simulations These days, the common opinion is that the flow and transport phenomena in a chemical reactor can be resolved in full detail across the whole domain under certain circumstances only. Laminar flow is an example of such a condition, as then the range of scales is limited. In addition, it is possible to describe or approximate the kinetics of the chemical reaction of interest by means of simplified rate expressions only. Van Wageningen et al.78 successfully simulated a waste-to-waste chemical process carried out in a static mixer and resulting in the formation of solid particles owing to fortunate conditions: the flow was unsteady and laminar, and the rate of the chemical reactions was transport-limited. By reducing the complex chemistry, Kleijn et al.55 were successful in reproducing deposition rates in chemical vapor deposition reactors of various types, which all are operated in the laminar flow regime. For processes carried out under turbulent-flow conditions, the spectrum of scales is broader, and simulating the fluid motions at all time and length scales concurrently still is virtually impossible. The limits of what is possible are shifting very fast these days, however. Recently, Chen et al.79 obtained an unprecedented, astonishing detailed picture of a gas-solid flow by carrying out a DNS by means of a macroscale particle method in which the gas phase was mimicked using 0.37 billion fluid particles, while the solid phase consisted of 0.1 million solid particles; the code was run on a platform comprising 256 graphical processing units (GPUs). Whether the flow has been resolved completely and is independent of grid size is unclear at the moment. For the time beingsand to avoid such excessive computational jobs as in the above paper due to Chen et al.79sone may think, however, of a multiscale approach combining a macroscale simulation in one way or another with one or more separate microscale simulations. This is certainly not a new idea: Ingram et al.45 already presented a framework for classifying different types of multiscale models. They illustrated their analysis by means of three different multiscale concepts for a catalytic packed bed reactor such as the partial oxidation of ortho-xylene to phthalic anhydride. They treated the reactor bulk fluid as the macroscopic domain (the reactor scale) and the catalyst pellets as the microscale domain. A uniform 1-D flow without axial dispersion and without radial temperature and concentration variations was assumed for the macroscopic domain. This latter assumption is a strong simplification, and actually the current paper exactly deals with the challenging task of resolving the detailed hydrodynamics both on the macroscale (the reactor) and on the microscale (e.g., the pellet scale). In this sense, the present paper is a follow-up of, e.g., the paper due to Lerou and Ng36 who made a plea for incorporating CFD in chemical reaction engineering modeling. Radl et al.89 also picked up this idea of incorporating the effect of hydrodynamics

Figure 6. Sample result from Radl et al.89 (their Figure 13c) showing the spatial distribution of the hydroxylamine concentration (nondimensional) around a large bubble cluster as calculated in a DNS for a periodic box for Re ) 38, Sc ) 50, Da1/Da2 ) 40, Ha ) 0.164, t* ) 81. Reprinted with permission from ref 89. Copyright 2008. Elsevier.

on chemical reactions: they studied fast chemical multiphase reactions in a mesoscale DNS of bubble swarms in a hydrogenation process. An impression of their results is presented in Figure 6. The latter paper followed an earlier paper by Khinast90 on fast chemical reactions in the wake of a single bubble. The present paper focuses on multiscale approaches in which the flow phenomena (fluid mechanics) are key and are resolved by means of the classical NS equation. Due to the strong nonlinearity of this equation, dedicated numerical solvers are essential. In the macroscale simulation, mesoscale or microscale effects on the full-scale processes (flow, turbulence, transport of heat and species, chemical reactions) are overlooked and modeled with a view of keeping the computational task within certain limits of size and time at the expense of accuracy. This is essentially what for single-phase flow systems has been done and is being done in the communities involved in modeling, e.g., nonreacting turbulent single-phase and two-phase flows,10,29,91,92 CVD reactors,53,55 and turbulent combustion.52,93 For fluidized beds, Kuipers has developed an extensive hierarchy of models (see, e.g., Van der Hoef et al.46-48) comprising several types of discrete particle and discrete bubble simulations (viz., resolved and unresolved types) and E-E simulations based on two-fluid models. The various types of computer simulations together form a family of multiscale simulations; the various simulations are carried out separately or at most sequentially. The discrete particle simulations are

Ind. Eng. Chem. Res., Vol. 49, No. 21, 2010

based on Newton’s law for colliding individual particles, with the interstitial gas phase being ignored and energy dissipation during the collisions being taken into account by means of a restitution coefficient and a friction coefficient representing the degree of inelasticity of the collisions.25 Such discrete particle simulations are essentially different from a Lagrangian approach in which particles respond to the forces exerted by the flowing carrier phase. While the latter approach applies to dilute particulate systems such as in gas-solid cyclones and solids suspensions in stirred vessels, dense fluidized beds may be dominated by particle collisions indeed. A rather simple way of incorporating the effect of microscale phenomena on the full-scale simulation is by means of phenomenological coefficients or constitutive models derived from off-line mesoscale simulations. Essentially, the idea is that separate mesoscale simulations are used to produce data and deduce closure laws to feed coarse-grained models which then can be used to compute the flow structures on the macroscale. Such mesoscale simulations are necessarily of the DNS type and are the subject of the next section. Kandhai et al.94 demonstrated the power of this approach by deriving the functional dependence of the single-particle drag force in a swarm of particles on the volume fraction by means of DNS of the fluid flow through disordered (or random) arrays of spheres in a periodic box. On the basis of these DNSs, these authors were able to assess for which values of the void fraction either the Wen and Yu correlation or the Ergun relation was most appropriate for modeling the fluid-particle drag in a coarsegrained simulation. Van der Hoef et al.47,95 and Beetstra et al.96 did something similar for both monodisperse and bidisperse particle assemblies. Sundaresan and co-workers used essentially the same procedure for deriving relations for specific fluidparticle interaction forces such as the lift force,97 the apparent virtual mass force,98 and the drag force.99 DNS in Periodic Boxes Let us now discuss the technique of carrying out a DNS for a small periodic box somewhere in the bulk of a reactor. The idea is that such a DNS may elucidate characteristic features of the processes of interest and that the increased understanding of the processes on the mesoscale may already result in more precise and more reliable predictions as to the outcome of the process of interest. The interest could be in the fluid-particle interaction force (see the above paragraphs) or in mass transport and chemical reactions in a swarm of bubbles. Another example is the mutual hydrodynamic interaction of a limited number of bubbles, drops, or particles including their readiness for collisions, break-up, and/or coalescence. Examples of such detailed studies by means of DNS are found in Ten Cate et al.,34 Sundaresan and co-workers,97-99 Tryggvason et al.,100 Derksen and Van den Akker,101,102 and Radl et al.89,103 Various numerical techniques are exploited in such studies, such as marker-and-cell, volume-of-fluid, front-tracking/ front-capturing, and lattice-Boltzmann. The reader is once more referred to Figure 6, from the work of Radl et al.89 who carried out a DNS for a periodic box containing a large bubble swarm. The presence of such a small box immersed in the dynamic ambiance of a turbulent flow is mimicked by using periodic boundary conditions, which protect the inside of the box against restricting effects of solid boundaries. By imposing typical (averaged) conditions of, e.g., flow, rate of energy dissipation, temperature, species concentrations, and/or volume fractions, we then are capable of studying how the flow and the various processes of interest evolve as a result of the governing models

10789

Figure 7. Cross-sectional area through a 3-D periodic box showing the results of a DNS of an assembly of colliding particles subjected to forced turbulence as reported by Ten Cate et al.34 (dots denote equally sized particles as cut through the cross-sectional plane. Not all velocity vectors are shown; colors are indicative of spatial variation in rate of energy dissipation). Reprinted with permission from ref 34. Copyright 2004. Cambridge University Press.

and equations. Whereas experimental techniques often fail in elucidating the mechanisms behind certain phenomena and processes, computational simulations are perfectly capable of doing this and hence a viable route for elucidating processes on the microscale level.79,89 This is a most welcome utilization of computational simulations indeed. Ten Cate et al.34 were interested in finding out locally about the frequencies of particle collisions in a suspension under the action of the turbulence of the carrier liquid. To this end, they performed a DNS of a particle suspension in a periodic box and used a specific forcing technique for subjecting the particles to turbulent flow conditions.104 The flow field around and between the interacting and colliding particles was fully resolved, while the particles were allowed to rotate in response to the vorticity in the surrounding carrier phase. Figure 7 illustrates how the flow field between colliding monosized spherical particles was completely resolved. Ten Cate et al.34 were able to learn from their DNSs about the mutual effect of microscale (particle scale) events and phenomena on the macroscale: the particle collisions are brought about by the turbulence, but the particles also modulate the turbulence. Energy spectra confirmed that the particles generate fluid motion at length scales on the order of the particle size. This results in a strong increase in the rate of energy dissipation on these length scales and in a decrease of the turbulence on larger length scales. A similar finding for an agitated liquid-liquid system was obtained by Derksen and Van den Akker.102 The simulations due to Ten Cate et al.34 may become more realistic when the method is extended to include particles of different sizes; it may be required to presume a certain particle size distribution thatsgiven the periodic boundariesshas to be retained during the DNS. In addition, it would be great when the method of running a DNS in a periodic box fed with local turbulence characteristics from a coarse-grained simulation could also be applied to fluid particles and would allow for bubble/

10790

Ind. Eng. Chem. Res., Vol. 49, No. 21, 2010

droplet coalescence and break-up in response to the local inbox turbulence level. Lattice Boltzmann Techniques In a remarkable development,105 lattice-Boltzmann (LB) techniques seem to offer a very attractive alternative to solving the classical NS equation for the purpose of tackling scale-bridging problems. The macroscopic continuum description embodied in the NS equation is replaced by a discrete particle-based, mesoscopicsor microscopic, though not truly molecularsview. Owing to their excellent numerical stability and constitutive versatility (see below), LB methods have started playing an important role as a simulation tool for advanced modeling of both materials properties and flow processes.58 Starting with Eggels106,107 and Derksen and co-workers,108,109 LB techniques are also being used for simulating macroscale turbulent flows. As a matter of fact, many of the above periodic-box DNSs were carried out by means of LB techniques as well. In the LB methodology,47,58,105,110 fluid flow is mimicked by a system of fictitious or pseudo-particles residing on a regular lattice and exhibiting discrete velocities which can be represented by some velocity distribution function. These pseudoparticles may be considered as coarse-grained groups of fluid molecules carrying some properties of real fluid portions. These pseudo-particles, or fluid portions, move at different speeds in specific, well-defined directions to neighboring lattice nodes. After having arrived on these nodes, the particles interact pairwise (“collide”) as hard spheres while particle number (mass), momentum, and energy are conserved; then they redistribute over the available velocity directions, which vary between the different LB versions. The steps of particles moving from node to node (translation, or streaming) and of particles colliding at a node (reaction) alternate. Several lattice topologies are in use, the face-centered hypercube (denoted as D3Q19) being one widely preferred.47,58,111 For the collision step, the Bhatnagar-Gross-Krook (BGK) model112 is quite common: it simplifies the complicated collisions to a relaxation process with a relaxation time constant, where it is assumed that the macroscopic quantities of interest are not influenced by the details of the collisions. With a single set of particles, a suitable lattice topology, and proper collision rules, and after a straightforward translation toward the fluid’s velocity and pressure fields, the single-phase fluid eventually obeys the classical NS equation (in the incompressible limit). Within a LB code, all information transfer is very local in time and in space (only to nearest neighbors). In this respect, LB essentially differs from the finite volume technique commonly used for solving the soft-field NS equation. As a result, LB codes are very suitable for massively parallel computations. Further, the LB method is superior due to the low number of operations per grid node and per time step. LB therefore classifies as a perfect technique for resolving by means of a DNS or a LES a broad spectrum of eddies in a turbulent flow as well as for multiscale strategies in general. One drawback of LB is that in principle it requires a homogeneous (regular) lattice of cubic cells.105 This is essentially different from the welcome feature of both finite element and finite volume methods, which exploit the idea of local mesh refinement in order to resolve strong local gradients, particularly in the vicinity of walls. In most LB applications, certainly for unsteady laminar and turbulent flows, however, fine meshes are intentionally used to resolve a substantial part of the spectrum of eddies. In addition, the inefficiency of using a uniform mesh is more than compensated by the high efficiency

per grid node. A particular though cumbersome type of grid refinement is achievable by combining several cubic and homogeneous meshes of different mesh sizes; of course, this requires matching of the different solutions at the corresponding interfaces by means of interpolation.61,113 Another aspect deserving special attention is the treatment of solid walls (stationary vessel walls, moving particles, revolving impeller) and the pertinent viscous wall layers. Rather than imposing a solid wall boundary condition, in many cases, an immersed boundary technique is to be preferred, just like in other computational techniques. In the particular immersed boundary technique used in our Delft group, the no-slip surfaces are represented by a set of points which do not need to coincide with the computational mesh. On each of these points, the noslip condition is enforced by adding momentum sources to the nearby grid points. The resulting body force field is confined to a layer of three grid spacings around the solid boundaries and is defined such that it drives the local fluid velocity toward the no-slip velocity at the wall. Further details of this immersed boundary method can be found elsewhere.108,114 Since LB deals with the properties of the fluid at the microkinetic level, it is perfectly capable of handling various thermodynamic properties of a fluid system. The latter is very handsome in view of multifluid systems. Not surprisingly, LB methods are exploited for modeling free surface flows and foaming.58,111 In the worldwide LB community, many are exploring and developing various LB techniques for simulating multiphase flows. A complete review of all this work is beyond the scope of the present paper. Of course, it is essential that the grid spacing of the computational grid for resolving the flow between the particles in a particle swarm is 1 order of magnitude smaller than the particle size. Some immersed boundary technique may enforce the zero-slip condition at the surface of the particles. For more details on simulating the multiphase flow of gas-solid suspensions, the reader is referred to some Delft papers.34,61,115 We will now focus on simulating liquid-liquid and gas-liquid dispersions where the particles are deformable, and coalescence and break-up are important and challenging issues. Lattice Boltzmann for Liquid-Liquid and Gas-Liquid Flows In the context of studying liquid-liquid dispersions in stirred vessels, Derksen and Van den Akker101,102 opted for a procedure due to He et al.116-118 based on a pressure distribution function, one for each of the two liquids which are represented by two sets of LB particles (say, red and blue). The method further uses an order parameter φ that stands for the local difference in the number densities of these red and blue particles. This φ tracks the interface between the two fluids; a transport equation for φ is solved by a LB method. The immiscibility of the two fluids is modeled by defining a free energy function that depends on φ and on the spatial gradients of φ; in this way, the effect of surface tension is mimicked. When properly chosen, this free energy function exhibits two minima which correspond with an equilibrium state preferred by fluid 1, and a second equilibrium state preferred by fluid 2. Derksen and Van den Akker101,102 were able to reproduce droplet break-up and coalescence within a two-fluid package in response to the typical varying turbulence levels such a package may encounter in a typical stirred vessel. Figure 8 illustrates the tendency toward coalescence when turbulence levels drop. In their simulations, the two liquids had the same density and viscosity however. The results of their exploratory

Ind. Eng. Chem. Res., Vol. 49, No. 21, 2010

10791

Figure 8. Snapshots from a movie due to Derksen and Van den Akker101 reproducing the evolution of the droplet size distribution in a periodic box as a result of an imposed decrease of turbulence level over time. Red color is indicative of droplets of liquid 1; blue denotes the immiscible carrier liquid. Liquids have the same density and viscosity.

study were qualitative and were not validated against quantitative experimental data, although reference is made to experiments reported by Schulze et al.119 Another aspect that requires further study relates to the thickness of the liquid-liquid interface that in this study extended over two or three lattice spacings, which was comparable with the Kolmogorov length scale. In physical reality, however, the interface is much thinner than any other length scale in turbulent liquid-liquid dispersions. Yet, their results looked very realistic! (This, however, is not uncommon, as many CFD approaches are successfully used beyond the bounds of their original validity range.52) As an alternative LB method for simulating multiphase flow, Shan and Chen120 propose to include in the LB formulation an intermolecular force, expressed as the gradient of a potential which is a function of the local fluid density. The form of the potential function dictates the resulting (nonideal) equation of state. When operated in the subcritical temperature range, the nonideal fluid comprises two coexisting states (phases) each having a different density. In addition to this density bifurcation, the interparticle or interaction potential provides a means of controlling the interfacial tension between the different phases by adding an attractive or repulsive tail to the elastic collisions of the fictitious particles making up the two phases. The concept of this interaction potential can be exploited to mimic drop formation and drop size in a gas-liquid mixture. Figure 9, from the work of Kamali et al.,121 shows how, starting from a random spatial distribution of gas and liquid “particles” (not shown), introducing the action of interaction potentials into the LB scheme results in a gradual phase separation and pertinent droplet formation. In the Shan and Chen LB method,120 however, the ratio of the densities of the two phases should be kept low as to avoid numerical instabilities (in Figure 9, this density ratio is still as low as 10.8). Kamali et al.121 succeeded in extending the Shan and Chen method toward higher density ratios, representative of a realistic gas-liquid system, by incorporating a different, more representative equation of state (EoS) than the ideal gas law implemented by Shan and Chen.120 In fact, introducing a different EoS exploits a very characteristic feature of LB, viz., its capability of dealing with intermolecular forces. This remedy against the low density ratio issue typical of the original Shan and Chen method was first suggested by Yuan and Schaefer.122 Following this suggestion, Yu et al.123 were able to get the density ratio as high as 60. Kamali et al.121 reported that introducing the Redlich-Kwong EoS resulted in a maximum density ratio of about 150, while the EoS due to Carnahan and Starling made a density ratio just beyond 1000 possible. Figure 10 is a snapshot from a movie obtained by Kamali et al.121 and illustrates the potential of the LB methods for reproducing the rise of a single spherical-cap-shaped gas bubble through a stagnant liquid inside a vertical tube. This result is

obtained by incorporating a hydrodynamic force on top of the interaction potential. So far, their results have not yet been validated against quantitative data. It would be very interesting and challenging to investigate the motion of multiple bubbles or drops responding to spatially and/or temporally varying turbulence levels in a carrier phase by adding proper forces into the LB recipe. A Concurrent Multiscale Strategy Pursuing the successes and promises of individual periodicbox DNSs which resolve local processes in all detail, we now propose to run a coarse-grained simulation and a number of DNSs concurrently. At regular time intervals, results from the two types of simulations are exchanged mutually. After each exchange of intermediate results, both types of simulations are continued until each of both types of simulations has converged. The essential step forward is that the results of a local DNS are fed back into the coarse-grained E-E two-phase flow simulation (RANS-based or LES) of a whole vessel or reactor. By this two-way coupling, local effects are also incorporated to a much larger extent, the impact of modeling is reduced, and more parameters are allowed to vary spatially. One option is to feed a new local value for ε back into the coarse-grained simulation, where this new local value of ε may be the result of local effects being mimicked in the concurrent DNS for a particular position in the vessel with the pertinent flow conditions. As discussed before, e.g., presence, motions, and collisions of dispersed-phase particles may modulate the turbulence in the carrier phase. A two-way coupling is essential to incorporate this turbulence modulation. Another option is to adapt bubble or drop size at the position for which the periodic-box DNS has been carried out on the basis of the turbulence characteristics found in the earlier phase of the coarse-grained simulation. The effect is that the particle size used in calculating the interaction force Fi in the two NS equations is allowed to respond locally to a higher or lower value of ε or shear rate than used so far in the coarse-grained simulation. The idea is then to continue again the coarse-grained simulation on the basis of the adapted bubble or drop size. Chemical reactors may provide a third examplesalso relevant for single-phase reactorssillustrating the relevance of carrying out a coarse-grained simulation concurrently with several (tens of) DNSs with a view to achieving frequent mutual exchange of intermediate results. In turbulent flows, the Batchelor scale governing chemical reactionsssee eq 2sis certainly smaller than the grid spacing of the coarse-grained simulation. Rather than considering the chemical reactions in the coarse-grained simulation, they might better be embedded in the periodic box DNS only. For slow reactions, micromixing models may be needed, such as the mechanistic cylindrical stretching vortex

10792

Ind. Eng. Chem. Res., Vol. 49, No. 21, 2010

Figure 10. A single gas bubble rising through a liquid, according to a LB simulation carried out by Kamali et al.121 Density ratio amounts to 1.326; the equation of state used is Van der Waals.

Figure 9. Snapshots at moments in time t ) 100, 200, 1000 from a sample calculation carried out by Kamali et al.121 and demonstrating the effect of a specific G value on the rate of phase separation in a LB two-phase simulation of the Shan and Chen type. Colors: red is liquid 1; blue is liquid 2. Density ratio is 10.8 only.

model proposed by Bakker and Van den Akker124 for a singlephase reactor. An alternative could be to represent the spatial distribution of all chemical species by means of fictitious particles carrying an array of concentrations, and to track these

fictitious particles within the periodic box with a Lagrangian Monte Carlo solver.41,42 Since most chemical reactions are accompanied by thermal effects, density and viscosity in the coarse-grained simulation of the chemical reactor may need frequent updating. When particles are formed as the result of the chemical reaction(s), their presence and motion may again modulate turbulence. Conversely, changes in the local variables such as density, viscosity, and turbulence intensity should be communicated to the DNS of the periodic boxes where the chemical reactions take place, as these variables may affect the local conditions. In particular for two-phase flow reactors, however, the proposed concurrent multiscale simulation strategy may be relevant, especially when one of the reactants is supplied via the gas bubbles, while the reaction itself takes place in the liquid phase. The DNS of the periodic box should then comprise both the mass transfer across the bubble interfaces (as affected by the flow around the bubbles) and the chemical reactions in the surrounding liquid. A nice example of such a periodic box DNS with fast chemical reactions is found in a recent paper from Radl et al.,89 who studied the hydrogenation of a nitroarene in a periodic box with a regular, fixed, 2-D grid and containing a “large” bubble cluster at a bubble Reynolds number of 38 and a gas hold-up of almost 5% (see Figure 6). These authors used a hybrid front-tracking/front-capturing method for dealing with the liquid-bubble interfaces. We think LB methods are faster and better geared to mimic the complex interaction of (turbulent) two-phase flow, mass (and heat) transfer, and chemical reactions. Of course, an important issue concerns the (necessarily limited) number W of periodic boxes for which a DNS has to be carried out that may be required for taking the spatial variations in shear rate, in ε, or in bubble or drop size across the full vessel into account. Selecting a proper number of periodic boxes and a proper spatial distribution of these boxes over the full vessel may be critical and may depend on the local gradients in shear rate or in ε. A related issue is the spatial interpolation across the full vessel between the varying drop or

Ind. Eng. Chem. Res., Vol. 49, No. 21, 2010

bubble sizes found from the various DNSs; such an interpolation may be necessary to avoid too large steps in bubble or drop size. Note that the periodic boxes need not be small compared to the vessel. It is very important that they contain a sufficiently large number of bubbles or dropssseen from a statistical point of viewswhile the hydrodynamic conditions inside the periodic box are inherently homogeneous. The only and eventual goal of the DNS for a periodic box is to come up with, e.g., a new bubble or drop size in response to the imposed turbulence level. The regular mutual exchange of simulation results deserves some more attention. After a (relatively large) number of iterations in the coarse-grained simulation, the two-phase flow field may have adapted more or less to the imposed spatial bubble or drop size distribution. Full convergence may not be required. Then, it may be time to feed updated turbulence characteristics to the W periodic boxes by means of the forcing technique and rerun the DNS for all W periodic boxes. This process is continued until convergence down to a suitable level has been attained in the coarse-grained simulation for the whole vessel. Ideally, the coarse-grained simulation and the various periodic box DNSs should be run truly concurrently. An alternative would be to make use of a look-up table connecting, e.g., bubble or drop size to ε (and possibly other variables). Such a look-up table may be constructed prior to running the coarse-grained simulation, or some in situ adaptive tabulation (ISAT) technique such as that exploited in simulating turbulent combustion125 might be useful here as well. The concurrent multiscale strategy for simulating two-phase flows advocated in this paper bears some resemblance to the hybrid formulations in use in the turbulent combustion field.52 Such hybrid formulations involve a LES combined with a lowdimensional (usually 1-D) model (such as a linear-eddy model or a 1-D turbulence model) for the subgrid scale processes. Such a low-dimensional model is there applied within each LES cell. The computational savings is attained by the reduced dimension. In the multiscale strategy proposed here for two-phase flow, the computational savings is realized by restricting the number W of periodic boxes for which a full 3-D DNS is carried out. As a matter of fact, the model advocated here is a typical example of dealing with what Echekki52 names “type-C coupling”: a strong physical coupling between processes while the scales are separated, requiring two-way coupling. The strategy proposed here is also an example of a concurrent multiscale approach as described by Ingram et al.,45 though on a by far more sophisticated level. The above concept of combining in a concurrent mode a coarse-grained E-E two-phase flow simulation for a whole vessel with DNSs for a number of periodic boxes is a truly multiscale approach. It is too optimistic to think that this concept will immediately result in a robust code that can be run routinely. In the first instance, finding out how exactly to combine the E-E simulation with several DNSs in a concurrent mode may be a tedious process. All such aspects as the frequency of the updating step, the number W of periodic boxes for which a DNS is carried out, their spatial distribution over the full vessel domain, and the interpolation in the full domain between the values fed back for, e.g., bubble or drop size, need to be sorted out. As part of the learning curve of how to run such a truly multiscale simulation, making use of a look-up table might be helpful. Combining a RANS-based simulation or preferably a LES for the whole vessel or reactor domain with several (ten or some

10793

tens) periodic-box DNSs in an iterative procedure really requires massive computer power. There is no a priori reason, however, whyscompared with the common approach of a constant bubble or drop sizesconvergence of the coarse-grained simulation would be more difficult in the newly proposed strategy. The new strategy aims at accommodating a spatial bubble or drop distribution which is more physical and realistic in cases where shear rate and ε vary spatially. Computational Issues Substantial progress has been made in exploiting CFD not just for simulating the turbulent flow field but also with the view of representing the various processes carried out in process equipment. In the early days of CFD, though, doubts were raised as to whether CFD would ever become a suitable tool for chemical engineers. In the meantime, however, due to the enormous growth in computer power and the proliferation of models and numerical methods, CFD has grown very powerful and versatile indeed. First of all, the increased computer power makes it possible to switch from steady-state simulations to transient simulations and to increase spatial resolution. One no longer has to be content with steady-flow RANS simulations. The number of grid points used for CFD studies of stirred tanks strongly varies with the type of CFD and generally increases with time. While for Bakker and Van den Akker91,92 the maximum number of grid points amounted to some 25 000, these days, hundreds of thousands of grid points are not uncommon for RANS simulations. Usually, LESs are carried out with the view of arriving at very fine spatial and temporal resolutions; then, substantially more grid points are desired than with RANS-based simulations. Full-scale LESs on fine grids of over 106 nodes currently are quite feasible and deliver realistic reproductions of transient flow and transport phenomena. Ten Cate et al.126 went as high as 35.5 × 106 grid points for their LES of an industrial crystallizer. Gillissen and Van den Akker66 carried out a DNS of a stirred vessel with baffles at Re ) 7300 on even a grid of some 2.9 × 109 grid cells. One may say that the size of the computational simulations keeps more or less equal pace with Moore’s law (saying the number of transistors on integrated circuits doubles every two years). There is no reason why computer power and CFDsin a type of leapfrog processswould not keep growing at the same pace. This means that what is advanced and time-consuming today will be a routine job in three or five years. Of course, academia may have the lead in this development. Furthermore, it has become attractive to speed up CFD simulations by running them on parallel computer platforms: a cluster of PCs operating under, e.g., LINUX or a massive parallel machine (e.g., a CRAY). In this way, larger CFD simulations can be run overnight or over the weekend. In addition, computer hardware (processor, memory) keeps becoming faster and cheaper. Vlachos et al.40 remarked that the low cost of Beowulf clusters renders multiscale simulation a reality. Just to illustrate the tempestuous growth in computational simulations: Derksen and Van den Akker109 carried out their first LESs of a stirred vessel on a four-processor HP machine where one impeller revolution took 11 h of wall-clock time. Derksen127 carried out his LES job on a Beowulf platform of 12 CPUs and needed only 7 h of wall-clock time for a single impeller revolution. In the huge (2.9 × 109 grid cells) DNS carried out by Gillissen and Van den Akker66 on 128 CPUs, one impeller revolution took over 33 h of wall-clock time due

10794

Ind. Eng. Chem. Res., Vol. 49, No. 21, 2010

to the very fine spatial resolution. As a matter of fact, such simulations yield incredible degrees of detail. Recently, graphical processing units (GPUs) have started being used for computational simulations to an increasing extent, as they turn out to be very fast. We are entering the era of supercomputing with Petaflops performance! Recently, one of the simulations of Chen et al.79 on a platform comprising 256 GPUs was 2600 times faster than on a single CPU core. The sky seems to be the limit. Last but not least, the LB methodology58,110 exhibits important advantages over conventional computational methods such as the common finite volume (FV) techniques. LB performs better in terms of computational speed due to the fact that all the information transfer is local in time and space: almost all communication is between nearest neighbor lattice nodes. As a result, LB is most suitable for massively parallel computing as it exhibits excellent scaling performance. Van Wageningen et al.62 and Van den Akker10 present comparisons between LB and FV (Fluent) simulations in terms of computational speed, which convincingly prove LB to be faster by at least 1 order of magnitude. Due to the requirement of robustness, a commercial CFD code should be capable of dealing with a wide variety of conditions; as a consequence, a commercial code is always in a disadvantageous position with respect to a dedicated research code developed in house. Yet, LB is inherently and evidently faster than a FV-based code. That is why LB is to be preferred for advanced, detailed, and time-consuming computational simulations: for DNSs of laminar single-phase flows in complex geometries, for LESs of turbulent single-phase flows, and for DNSs of multiphase flows. Conclusions This review has reported on the wealth of computational fluid dynamics approaches and computational multiscale models developed in, say, the past 30 years. Substantial progress has been realized, not in the least due to the tempestuous growth in computational resources and the development of very accurate and powerful numerical algorithms. Grid size, job size, computational time, processor speed, computer memory: all have increased substantially. As a result, the temporal and spatial resolutions keep increasing and are delivering ever more unprecedented amounts of detail as to the local microscale or mesoscale flow and transport phenomena. An inventory has been made of quite a few multiscale modeling approaches developed in such divergent fields as chemical reactor engineering, turbulence modeling, combustion, chemical vapor deposition, and materials science. Although the challenges are identical and in some way the remedies also bear a strong resemblance, it strongly depends on the kind of phenomena and processes at stake as to which methods or equations are used for describing them. Models describing single-phase and two-phase flow involve the Navier-Stokes equation, and this turns processes involving flows differentlysat least at a first glance. The essentials of direct numerical simulations (DNSs), RANS-based simulations, large eddy simulations (LESs), Euler-Lagrange simulations, and Euler-Euler (or two-fluid) simulations have been summarized. While mostly finite volume techniques are used for solving the highly nonlinear NS equation, which is at the heart of all these flow simulations, in the past decade, lattice Boltzmann (LB) techniques have proven to be equally or even better suited. This might reduce the difference between the multiscale models used in the various

fields: they all may be described by means of particle-based simulation techniques indeed. The focus of this paper is on the results and promises of carrying out DNSs for the processes on the small scale at some position in the whole-flow domain. By means of a suitable forcing technique, the turbulent flow conditions as prevailing at the position of interest are fed to a box with periodic walls; the DNS then resolves the flow in all detail. The use of periodic boundary conditions aims at representing a box immersed in and subjected to the dynamic ambiance of the turbulent carrier phase. The periodic box may contain particles, drops, or bubbles which may interact directly (during collisions) or via the carrier phase. Such a DNS provides a forum for two-way coupling: the carrier phase affects particle motion, and the particles affect the fluid’s motion. LB techniques have proven to be very suitable and fast for such simulations. Throughout the paper, several promising results have been presented from LB simulations of the DNS type for two-phase flows. Different LB models for dealing with two fluid phases and the phase interface are available and have been implemented in our in-house code. An important issue is the ratio of the two fluid densities since most LB methods so far suffered from the restriction of yielding converged results for low density ratios only. Much is to be awaited from LB methods provided that effective expressions are found for dealing with the interaction between the two phases and for the equation of state of the individual fluid phases. The idea is that properly designed LB methods will soon become capable of reproducing the effect of turbulence on bubble and drop size (distributions). In the literature, quite a few papers can be identified which all implicitly or explicitly aim at reducing the role of closure models representing unresolved physics or chemistry. Several papers in the literature report direct numerical simulations either for the full domain of a reactor or for a periodic box somewhere in the domain. These simulations throw new light on the intimate processes in turbulent single-phase and two-phase flows. Some authors try and derive from their direct numerical simulations simple correlations which can be used in coarse-grained simulations for the whole reactor. Others have built a hierarchy of models which are run separately or sequentially. A truly multiscale approach has been proposed for a gas-liquid or liquid-liquid system, combining in a concurrent mode a coarse-grained Euler-Euler two-phase flow simulation (RANS-based or LES) for a whole vessel with DNSs for a number of periodic boxes. The idea is that the coarse-grained simulation of the whole vessel receives at regular time intervals feedback about local bubble or drop sizes or even their distribution. This local bubble or drop size information should result from running DNSs for several periodic boxes subject to the different turbulent flow characteristics pertinent to specific positions in the vessel. In an iterative procedure that makes use of a frequent exchange of intermediate simulation results between the two simulation types, convergence should be reached toward a coarse-grained simulation result exhibiting a spatial bubble or drop size distribution. The concept of concurrently running a coarse-grained simulation and several (some tens?) periodic-box DNSssincluding periodic two-way exchange of intermediate results until convergence has been reachedsalso attempts to reduce the role of closure models but in quite a different way than before. The concept differs from the common Euler-Euler approach in which a single bubble or drop size is prescribed for the whole vessel. It also differs from Euler-Euler simulations extended with methods of moments representing particle size distributions.

Ind. Eng. Chem. Res., Vol. 49, No. 21, 2010

Further, there is the option to include chemical reactions in the periodic-box DNS. At the moment, this concept of concurrent multiscale simulations has not yet been implemented, but several issues and challenges have been addressed in this paper. Literature Cited (1) Levenspiel, O. Chemical Reaction Engineering, 3rd ed.; Wiley: New York, 1999. (2) Villermaux, J. Genie de la reaction chimique, 2nd revised ed.; Lavoisier: Paris, 1995. (3) Godbole, S. P., Shah, Y. T. Design and operation of bubble column reactors. In Encyclopedia of Fluid Mechanics: Vol. 3 - Gas-Liquid Flows; Cheremisinoff, N. P., Ed.; Gulf Publishing Company: Houston, TX, 1986. (4) Mills, P. L.; Chaudhari, R. V. Multiphase catalytic reactor engineering and design for pharmaceuticals and fine chemicals. Catal. Today 1997, 37 (4), 367–404. (5) Van Geem, K. M.; Heynderickx, G. J.; Marin, G. B. Effect of Radial Temperature Profiles on Yields in Steam Cracking. AIChE J. 2004, 50 (1), 173–183. (6) Van Deemter, J. J. Mixing and Contacting in Gas Solid Fluidized Beds. Chem. Eng. Sci. 1961, 13 (3), 143–154. (7) Kuipers, J. A. M.; Van Swaaij, W. P. M. Application of computational fluid dynamics to chemical reaction engineering. ReV. Chem. Eng. 1997, 13 (3), 1–118. (8) Van den Akker, H. E. A., Computational Fluid Dynamics: More Than a Promise To Chemical Reaction Engineering. In CHISA, Keynote Paper #1270, Prague, CZ, 2000. (9) Van den Akker, H. E. A. Keynote speech: On status and merits of computational fluid dynamics. In 4th International Conference on Bioreactor and Bioprocess Fluid Dynamics; BHR Group: Edinburgh, U.K., 1997; Vol. 25, pp 407-432. (10) Van den Akker, H. E. A. The Details of Turbulent Mixing Processes and their Simulation. In AdVances in Chemical Engineering; Elsevier Academic Press: London, 2006; Vol 31, pp 151-229. (11) Chen, R. C.; Reese, J.; Fan, L. S. Flow structure in a 3-dimensional bubble column and 3-phase fluidized bed. AIChE J. 1994, 40 (7), 1093– 1104. (12) Delnoij, E.; Kuipers, J. A. M.; Van Swaaij, W. P. M.; Westerweel, J. Measurement of gas-liquid two-phase flow in bubble columns using ensemble correlation PIV. Chem. Eng. Sci. 2000, 55 (17), 3385–3395. (13) Mudde, R. F.; Groen, J. S.; Van den Akker, H. E. A. Liquid velocity field in a bubble column: LDA experiments. Chem. Eng. Sci. 1997, 52 (2122), 4217–4224. (14) Mudde, R. F.; Groen, J. S.; Van den Akker, H. E. A. Application of LDA to bubbly flows. Nucl. Eng. Des. 1998, 184 (2-3), 329–338. (15) Groen, J. S.; Mudde, R. F.; Van den Akker, H. E. A. Time dependent behaviour of the flow in a bubble column. Chem. Eng. Res. Des. 1995, 73 (A6), 615–621. (16) Groen, J. S.; Oldeman, R. G. C.; Mudde, R. F.; Van den Akker, H. E. A. Coherent structures and axial dispersion in bubble column reactors. Chem. Eng. Sci. 1996, 51 (10), 2511–2520. (17) Julia, J. E.; Harteveld, W. K.; Mudde, R. F.; Van den Akker, H. E. A. On the accuracy of the void fraction measurements using optical probes in bubbly flows. ReV. Sci. Instrum. 2005, 76 (3), 1–13. (18) Xue, J.; Al-Dahhan, M.; Dudukovic, M. P.; Mudde, R. F. Bubble velocity, size, and interfacial area measurements in a bubble column by four-point optical probe. AIChE J. 2008, 54 (2), 350–363. (19) Harteveld, W. K.; Mudde, R. F.; Van den Akker, H. E. A. Dynamics of a bubble column: Influence of gas distribution on coherent structures. Can. J. Chem. Eng. 2003, 81 (3-4), 389–394. (20) Li, J. H.; Kwauk, M. Multiscale nature of complex fluid-particle systems. Ind. Eng. Chem. Res. 2001, 40 (20), 4227–4237. (21) Mudde, R. F. Gravity-driven bubbly flows. Annu. ReV. Fluid Mech. 2005, 37, 393–423. (22) Van den Akker, H. E. A. Coherent structures in multiphase flows. Powder Technol. 1998, 100 (2-3), 123–136. (23) Mudde, R. F.; Saito, T. Hydrodynamical similarities between bubble column and bubbly pipe flow. J. Fluid Mech. 2001, 437, 203–228. (24) Bauer, M.; Eigenberger, G. A concept for multi-scale modeling of bubble columns and loop reactors. Chem. Eng. Sci. 1999, 54 (21), 5109– 5117. (25) Hoomans, B. P. B.; Kuipers, J. A. M.; Briels, W. J.; Van Swaaij, W. P. M. Discrete particle simulation of bubble and slug formation in a two-dimensional gas-fluidised bed: A hard-sphere approach. Chem. Eng. Sci. 1996, 51 (1), 99–118. (26) Verloop, W. C. B. D.; Van den Akker, H. E. A.; Hein, K. R. G. The fluid dynamics of particles in the freeboard of a pressurized fluidized

10795

bed combustor. In 12th International Conference on Fluidized Bed Combustion, San Diego, CA; Berkely Electronic Press: Berkely, CA, 1993; pp 53-62. (27) Verloop, W. C.; Van den Akker, H. E. A.; Boersma, D.; Hein, K. R. G. Transient modelling of the combined gas-particle flow in the freeboard of a fluidized bed reactor. Proceedings of the International Conference on Fluidized Bed Combustion; Berkely Electronic Press: Berkely, CA, 1995; pp 789-799. (28) Mudde, R. F.; Harteveld, W. K.; Van den Akker, H. E. A.; Van Der Hagen, T. H. J. J.; Van Dam, H. Gamma radiation densitometry for studying the dynamics of fluidized beds. Chem. Eng. Sci. 1999, 54 (1314), 2047–2054. (29) Agrawal, K.; Loezos, P. N.; Syamlal, M.; Sundaresan, S. The role of meso-scale structures in rapid gas-solid flows. J. Fluid Mech. 2001, 445, 151–185. (30) Tennekes, H.; Lumley, J. L. A First Course in Turbulence; MIT Press: Cambridge, MA, 1972. (31) Van den Akker, H. E. A. The Euler-Euler approach to dispersed two-phase flows in the turbulent regime. ERCOFTAC Bull. 1998, 36, 30– 33. (32) Crowe, C. T. On models for turbulence modulation in fluid-particle flows. Int. J. Multiphase Flow 2000, 26, 719–727. (33) Sato, Y.; Fukuichi, U.; Hishida, K. Effect of inter-particle spacing on turbulence modulation by Lagrangian PIV. Int. J. Heat Fluid Flow 2000, 21, 554–561. (34) Ten Cate, A.; Derksen, J. J.; Portela, L. M.; Van den Akker, H. E. A. Fully resolved simulations of colliding monodisperse spheres in forced isotropic turbulence. J. Fluid Mech. 2004, 519, 233–271. (35) Deen, W. M. Analysis of Transport Phenomena; Oxford University Press: New York, 1998. (36) Lerou, J. J.; Ng, K. M. Chemical reaction engineering: A multiscale approach to a multiobjective task. Chem. Eng. Sci. 1996, 51 (10), 1595– 1614. (37) Li, J.; Ge, W.; Zhang, J.; Kwauk, M. Multi-scale compromise and multi-level correlation in complex systems. Chem. Eng. Res. Des. 2005, 83 (A6), 574–582. (38) Wei, J. Contribution to General discussion on kinetics, catalysis, and reactor engineering. AdVances in Chemical Engineering; Colton, C. K., Ed.; Elsevier Academic Press: London, 1991; Vol 16, p 253. (39) Prasad, V.; Vlachos, D. G. Multiscale model and informatics-based optimal design of experiments: Application to the catalytic decomposition of ammonia on ruthenium. Ind. Eng. Chem. Res. 2008, 47 (17), 6555– 6567. (40) Vlachos, D. G.; Mhadeshwar, A. B.; Kaisare, N. S. Hierarchical multiscale model-based design of experiments, catalysts, and reactors for fuel processing. Comput. Chem. Eng. 2006, 30 (10-12), 1712–1724. (41) Van Vliet, E.; Derksen, J. J.; Van den Akker, H. E. A. Turbulent mixing in a tubular reactor: Assessment of an FDF/LES approach. AIChE J. 2005, 51 (3), 725–739. (42) Van Vliet, E.; Derksen, J. J.; Van den Akker, H. E. A.; Fox, R. O. Numerical study on the turbulent reacting flow in the vicinity of the injector of an LDPE tubular reactor. Chem. Eng. Sci. 2007, 62 (9), 2435–2444. (43) Petre, C. F.; Larachi, F.; Iliuta, I.; Grandjean, B. P. A. Pressure drop through structured packings: Breakdown into the contributing mechanisms by CFD modeling. Chem. Eng. Sci. 2003, 58 (1), 163–177. (44) Dixon, A. G.; Nijemeisland, M.; Stitt, E. H. Packed tubular reactor modeling and catalyst design using computational fluid dynamics. In AdVances in Chemical Engineering; Marin, G. B. , Ed.; Elsevier Academic Press: London, 2006; Vol. 31, pp 307-389. (45) Ingram, G. D.; Cameron, I. T.; Hangos, K. M. Classification and analysis of integrating frameworks in multiscale modelling. Chem. Eng. Sci. 2004, 59, 2171–2187. (46) Van der Hoef, M. A.; Van Sint Annaland, M.; Kuipers, J. A. M. Computational fluid dynamics for dense gas-solid fluidized beds: A multiscale modeling strategy. Chem. Eng. Sci. 2004, 59 (22-23), 5157–5165. (47) Van der Hoef, M. A.; Ye, M.; Van Sint Annaland, M.; Andrews, A. T.; Sundaresan, S.; Kuipers, J. A. M. Multiscale Modeling of GasFluidized Beds. In AdVances in Chemical Engineering; Elsevier Academic Press: London, 2006; Vol. 31, pp 65-149. (48) Van der Hoef, M. A.; Van Sint Annaland, M.; Deen, N. G.; Kuipers, J. A. M. Numerical simulation of dense gas-solid fluidized beds: A multiscale modeling strategy. In Annual ReView of Fluid Mechanics; Annual Reviews: Palo Alto, CA, 2008; Vol. 40, pp 47-70. (49) Wang, J. A review of Eulerian simulation of Geldart A particles in gas-fluidized beds. Ind. Eng. Chem. Res. 2009, 48 (12), 5567–5577. (50) Wang, W.; Lu, B.; Zhang, N.; Shi, Z. S.; Li, J. H. A review of multiscale CFD for gas-solid CFB modeling. Int. J. Multiphase Flow 2010, 36, 109–118.

10796

Ind. Eng. Chem. Res., Vol. 49, No. 21, 2010

(51) Peters, N. Multiscale combustion and turbulence. Proc. Combust. Inst. 2009, 32, 1–25. (52) Echekki, T. Topical review: Multiscale methods in turbulent combustion: strategies and computational challenges. Comput. Sci. DiscoVery 2009, 2, 013001. (53) Rodgers, S. T.; Jensen, K. F. Multiscale modeling of chemical vapor deposition. J. Appl. Phys. 1998, 83 (1), 524–530. (54) Kleijn, C. R. Computational modeling of transport phenomena and detailed chemistry in chemical vapor deposition - a benchmark solution. Thin Solid Films 2000, 365 (2), 294–306. (55) Kleijn, C. R.; Dorsman, R.; Kuijlaars, K. J.; Okkerse, M.; Van Santen, H. Multi-scale modeling of chemical vapor deposition processes for thin film technology. J. Cryst. Growth 2007, 303 (1 SPEC. ISS.), 362– 380. (56) Abbate, G.; Thijsse, B. J.; Kleijn, C. R. Validation of a hybrid Navier-Stokes/DSMC method for multiscale transient and steady-state gas flows. Int. J. Multiscale Comput. Eng. 2008, 6 (1), 1–12. (57) Chatterjee, A.; Vlachos, D. G. An overview of spatial microscopic and accelerated kinetic Monte Carlo methods. J. Comput.-Aided Mater. Des. 2007, 14 (2), 253–308. (58) Raabe, D. Overview of the lattice Boltzmann method for nanoand microscale fluid dynamics in materials science and engineering. Modell. Simul. Mater. Sci. Eng. 2004, 12 (6), R13–R46. (59) Kleijn, C. R., Numerical simulation of flow and chemistry in thermal chemical vapor deposition processes. In Proceedings of the NATO AdVanced Study Institute on Chemical Physics of Thin Film Deposition Processes for Micro- and Nanotechnologies,, Kaunas, Lithuania; Pauleau, Y., Ed.; Kluwer Academic Publishers: Dordrecht, The Netherlands, 2001; pp 119-144. (60) Ten Cate, A.; Nieuwstad, C. H.; Derksen, J. J.; Van den Akker, H. E. A. Particle imaging velocimetry experiments and lattice-Boltzmann simulations on a single sphere settling under gravity. Phys. Fluids 2002, 14 (11), 4012–4025. (61) Rohde, M.; Derksen, J. J.; Van den Akker, H. E. A. An applicability study of advanced lattice-Boltzmann techniques for moving, no-slip boundaries and local grid refinement. Comput. Fluids 2008, 37 (10), 1238– 1252. (62) Van Wageningen, W. F. C.; Kandhai, D.; Mudde, R. F.; Van den Akker, H. E. A. Dynamic flow in a Kenics static mixer: An assessment of various CFD methods. AIChE J. 2004, 50 (8), 1684–1696. (63) Wilcox, D. C. Turbulence Modelling for CFD; DCW Industries Inc.: La Canada, CA, 1993. (64) Moin, P.; Kim, J. Tackling turbulence with supercomputers. Sci. Am. 1997, 276 (1), 62–68. (65) Sommerfeld, M.; Decker, S. State of the art and future trends in CFD simulation of stirred vessel hydrodynamics. Chem. Eng. Technol. 2004, 27 (3), 215–224. (66) Gillissen, J. J. J.; Van den Akker, H. E. A. Direct numerical simulation of the turbulent flow in a baffled tank driven by a Rushton turbine. Presented at Mixing XXII (North American Mixing Forum), June 2010. (67) Frisch, U. Turbulence - The Legacy of A. N. KolmogoroV; Cambridge University Press: Cambridge, U.K., 1995. (68) Rodi, W. Turbulence models and their application in hydraulics a state of the art reView; International Association for Hydraulic Research: Delft, The Netherlands, 1984. (69) Rielly, C. D.; Marquis, A. J. A particle’s eye view of crystallizer fluid mechanics. Chem. Eng. Sci. 2001, 56 (7), 2475–2493. (70) Smagorinsky, J. General circulation experiments with the primitive equations: 1. The basic experiment. Monthly Weather ReV. 1963, 91, 99– 164. (71) Derksen, J. Assessment of large eddy simulations for agitated flows. Chem. Eng. Res. Des. 2001, 79 (A8), 824–830. (72) Derksen, J. J. Numerical simulation of solids suspension in a stirred tank. AIChE J. 2003, 49 (11), 2700–2714. (73) Hartmann, H.; Derksen, J. J.; Montavon, C.; Pearson, J.; Hamill, I. S.; Van den Akker, H. E. A. Assessment of large eddy and RANS stirred tank simulations by means of LDA. Chem. Eng. Sci. 2004, 59 (12), 2419– 2432. (74) Hoekstra, A. J.; Derksen, J. J.; Van den Akker, H. E. A. An experimental and numerical study of turbulent swirling flow in gas cyclones. Chem. Eng. Sci. 1999, 54 (13-14), 2055–2065. (75) Derksen, J. J.; Van den Akker, H. E. A.; Sundaresan, S. Two-way coupled large-eddy simulations of the gas-solid flow in cyclone separators. AIChE J. 2008, 54 (4), 872–885. (76) Patankar, N. A.; Joseph, D. D. Modeling and numerical simulation of particulate flows by the Eulerian-Lagrangian approach. Int. J. Multiphase Flow 2001, 27 (10), 1659–1684. (77) Hartmann, H.; Derksen, J. J.; Van den Akker, H. E. A. Numerical simulation of a dissolution process in a stirred tank reactor. Chem. Eng. Sci. 2006, 61 (9), 3025–3032.

(78) Van Wageningen, W. F. C.; Mudde, R. F.; Van den Akker, H. E. A. Numerical simulation of growing Cu particles in a Kenics static mixer reactor in which Cu2+ is reduced by carbohydrates. Chem. Eng. Sci. 2004, 59 (22-23), 5193–5200. (79) Chen, F. G.; Ge, W.; He, X. F.; Li, B.; Li, J. H.; Li, X. P.; Wang, X. W.; Yuan, X. L. Multi-scale HPC system for multi-scale discrete simulation-Development and application of a supercomputer with 1 Petaflops peak performance in single precision. Particuology 2009, 7, 332–335. (80) Rietema, K.; Van den Akker, H. E. A. On the momentum equations in dispersed two-phase systems. Int. J. Multiphase Flow 1983, 9 (1), 21– 36. (81) Oey, R. S.; Mudde, R. F.; Van den Akker, H. E. A. Sensitivity study on interfacial closure laws in two-fluid bubbly flow simulations. AIChE J. 2003, 49 (7), 1621–1636. (82) Zhou, L. X. Advances in studies on two-phase turbulence in dispersed multiphase flows. Int. J. Multiphase Flow 2010, 36 (2), 100– 108. (83) Lance, M.; Marie, J. L.; Bataille, J. Homogeneous turbulence in bubbly flows. J. Fluids Eng.sTrans. ASME 1991, 113 (2), 295–300. (84) Curtis, J. S.; Van Wachem, B. Modeling particle-laden flows: A research outlook. AIChE J. 2004, 50 (11), 2638–2645. (85) Portela, L. M.; Oliemans, R. V. A. Possibilities and limitations of computer simulations of industrial turbulent dispersed multiphase flows. Flow Turbulence Combust. 2006, 77 (1-4), 381–403. (86) Niceno, B.; Dhotre, M. T.; Deen, N. G. One-equation sub-grid scale (SGS) modelling for Euler-Euler large eddy simulation (EELES) of dispersed bubbly flow. Chem. Eng. Sci. 2008, 63 (15), 3923–3931. (87) Venneker, B. C. H.; Derksen, J. J.; Van den Akker, H. E. A. Population balance modeling of aerated stirred vessels based on CFD. AIChE J. 2002, 48 (4), 673–685. (88) Petitti, M.; Nasuti, A.; Marchisio, D. L.; Vanni, M.; Baldi, G.; Mancini, N.; Podenzani, F. Bubble Size Distribution Modeling in Stirred Gas-Liquid Reactors with QMOM Augmented by a New Correction Algorithm. AIChE J. 2010, 56 (1), 36–53. (89) Radl, S.; Koynov, A.; Tryggvason, G.; Khinast, J. G. DNS-based prediction of the selectivity of fast multiphase reactions: Hydrogenation of nitroarenes. Chem. Eng. Sci. 2008, 63 (12), 3279–3291. (90) Khinast, J. G. Impact of 2-D bubble dynamics on the selectivity of fast gas-liquid reactions. AIChE J. 2001, 47, 2304–2319. (91) Bakker, A.; Van den Akker, H. E. A. Single-phase flow in stirred reactors. Chem. Eng. Res. Des. 1994, 72 (A4), 583–593. (92) Bakker, A.; Van den Akker, H. E. A. Computational model for the gas-liquid flow in stirred reactors. Chem. Eng. Res. Des. 1994, 72 (A4), 594–606. (93) Chen, J. C.; Choudhary, A.; De Supinski, B.; DeVries, M.; Hawkes, E. R.; Klasky, S.; Liao, W. K.; Ma, K. L.; Mellor-Crummey, J.; Podhorszki, N.; Sankaran, R.; Shende, S.; Yoo, C. S. Terascale direct numerical simulations of turbulent combustion using S3D. Comput. Sci. DiscoVery 2009, 2, 015001. (94) Kandhai, D.; Derksen, J. J.; Van den Akker, H. E. A. Interphase drag coefficients in gas-solid flows. AIChE J. 2003, 49 (4), 1060–1065. (95) Van der Hoef, M. A.; Beetstra, R.; Kuipers, J. A. M. LatticeBoltzmann simulations of low-Reynolds-number flow past mono- and bidisperse arrays of spheres: Results for the permeability and drag force. J. Fluid Mech. 2005, 528, 233–254. (96) Beetstra, R.; Van der Hoef, M. A.; Kuipers, J. A. M. Drag force of intermediate Reynolds number flow past mono- and bidisperse arrays of spheres. AIChE J. 2007, 53 (2), 489–501. (97) Sankaranarayanan, K.; Sundaresan, S. Lift force in bubbly suspensions. Chem. Eng. Sci. 2002, 57 (17), 3521–3542. (98) Ten Cate, A.; Sundaresan, S. Analysis of unsteady forces in ordered arrays of monodisperse spheres. J. Fluid Mech. 2006, 552, 257–287. (99) Holloway, W.; Yin, X. L.; Sundaresan, S. Fluid-Particle Drag in Inertial Polydisperse Gas-Solid Suspensions. AIChE J. 2010, 56 (8), 1995– 2004. (100) Tryggvason, G.; Esmaeeli, A.; Lu, J. C.; Biswas, S. Direct numerical simulations of gas/liquid multiphase flows. Fluid Dynamics Res. 2006, 38 (9), 660–681. (101) Derksen, J. J., Van den Akker, H. E. A. Multi-scale simulations of stirred liquid-liquid dispersions. In 12th European Conference on Mixing; Magelli, F., Baldi, G., Brucato, A., Eds.; AIDIC: Bologna, Italy, 2006 [CDROM]. (102) Derksen, J. J.; Van den Akker, H. E. A. Multi-scale simulations of stirred liquid-liquid dispersions. Chem. Eng. Res. Des. 2007, 85 (5 A), 697–702. (103) Radl, S.; Tryggvason, G.; Khinast, J. G. Flow and mass transfer of fully resolved bubbles in non-Newtonian fluids. AIChE J. 2007, 53 (7), 1861–1878.

Ind. Eng. Chem. Res., Vol. 49, No. 21, 2010 (104) Ten Cate, A.; Van Vliet, E.; Derksen, J. J.; Van den Akker, H. E. A. Application of spectral forcing in lattice-Boltzmann simulations of homogeneous turbulence. Comput. Fluids 2006, 35 (10), 1239–1251. (105) Succi, S. Lattice Boltzmann equation: failure or success. Physica A 1997, 240, 221–228. (106) Eggels, J. G. M. Direct and large-eddy simulation of turbulent fluid flow using the lattice-Boltzmann scheme. Int. J. Heat Fluid Flow 1996, 17 (3), 307–323. (107) Eggels, J. G. M. Large-eddy simulation of turbulent flows in baffled stirred tank reactors. In AdVances in Turbulences VI; Gavrilakis, S., Machiels, L., Monkewitz, P. A., Eds.; Infoscience: Lausanne, Switzerland, 1996; Vol. 36, pp 19-22. (108) Derksen, J. J.; Kooman, J. L.; Van den Akker, H. E. A. Parallel fluid flow simulations by means of a Lattice Boltzmann scheme. Lect. Notes Comput. Sci. 1997, 1225, 542. (109) Derksen, J.; Van den Akker, H. E. A. Large eddy simulations on the flow driven by a Rushton turbine. AIChE J. 1999, 45 (2), 209–221. (110) Succi, S. The Lattice Boltzmann Equation for Fluid Dynamics and Beyond; Oxford University Press: New York, 2001. (111) Korner, C.; Thies, M.; Hofmann, T.; Thurey, N.; Rude, N. Lattice Boltzmann model for free surface flow for modeling foaming. J. Stat. Phys. 2005, 121 (1/2), 179–196. (112) Bhatnagar, P. L.; Gross, E. P.; Krook, M. A model for collision processes in gases. Phys. ReV. 1954, 94, 511. (113) Rohde, M.; K, D.; Derksen, J. J.; Van den Akker, H. E. A. A generic, mass conservative local grid refinement technique for latticeBoltzmann schemes. Int. J. Numer. Meth. Fluids 2006, 51 (4), 439–468. (114) Goldstein, D.; Handler, R.; Sirovich, L. Modeling a No-Slip Flow Boundary with an External Force Field. J. Comput. Phys. 1993, 105, 354– 366. (115) Rohde, M.; Kandhai, D.; Derksen, J. J.; Van den Akker, H. E. A. Improved bounce-back methods for no-slip walls in lattice-Boltzmann schemes: Theory and simulations. Phys. ReV. E: Stat., Nonlinear, Soft Matter Phys. 2003, 67 (6 2), 066703/1–066703/10. (116) He, X. Y.; Zhang, R. Y.; Chen, S. Y.; Doolen, G. D. On the threedimensional Rayleigh-Taylor instability. Phys. Fluids 1999, 11 (5), 1143– 1152.

10797

(117) He, X. Y.; Chen, S. Y.; Zhang, R. Y. A lattice Boltzmann scheme for incompressible multiphase flow and its application in simulation of Rayleigh-Taylor instability. J. Comput. Phys. 1999, 152 (2), 642–663. (118) Zhang, R. Y.; He, X. Y.; Chen, S. Y. Interface and surface tension in incompressible lattice Boltzmann multiphase model. Comput. Phys. Commun. 2000, 129 (1-3), 121–130. (119) Schulze, K. R. J.; Kraume, M., Investigations of local drop size distributions and scale-up in stirred liquid-liquid dispersions. In 10th European Conference on Mixing; Van den Akker, H. E. A., Derksen, J. J., Eds.; Elsevier: Delft, The Netherlands, 2000; pp 181-188. (120) Shan, X. W.; Chen, H. D. Lattice Boltzmann model for simulating flows with multiple phases and components. Phys. ReV. E 1993, 47 (3), 1815–1819. (121) Kamali, M. R., Sundaresan, S., Gillissen, J. J. J., Van den Akker, H. E. A. Manuscript in preparation, 2010. (122) Yuan, P.; Schaefer, L. Equations of state in a lattice Boltzmann model. Phys. Fluids 2006, 18, 4. (123) Yu, Z.; Heraminger, O.; Fan, L. S. Experiment and lattice Boltzmann simulation of two-phase gas-liquid flows in microchannels. Chem. Eng. Sci. 2007, 62 (24), 7172–7183. (124) Bakker, R. A.; Van den Akker, H. E. A. A Lagrangian description of micromixing in a stirred tank reactor using 1D-micromixing models in a CFD flow field. Chem. Eng. Sci. 1996, 51 (11), 2643–2648. (125) Pope, S. B. PDF Methods for turbulent reactive flows. Prog. Energy Combust. Sci. 1985, 11 (2), 119–192. (126) Ten Cate, A.; Bermingham, S. K.; Derksen, J. J.; Kramer, H. M. J. Compartmental modeling of an 1100L DTB crystallizer based on large eddy flow simulation. In 10th European Conference on Mixing; Van den Akker, H. E. A., Derksen, J. J., Eds.; Elsevier: New York, 2000; pp 255-264. (127) Derksen, J. J. Long-time solids suspension simulations by means of a large-eddy approach. Chem. Eng. Res. Des. 2006, 84 (A1), 38–46.

ReceiVed for reView March 15, 2010 ReVised manuscript receiVed August 9, 2010 Accepted August 25, 2010 IE1006382