Pergamon Chemical Engineering Science Harris Amsterdam

Pergamon Chemical Engineering Science, Vol. 51, No. 10, pp. 1569-1594, 1996 Copyright © 1996 Elsevier Science Ltd Printed in Great Britain. All rights reserved S0009-2509(96)00021-8 0009-2509/96 $15.00 + 0.130 Computational Fluid Dynamics for Chemical Reactor Engineering C.K. Harris, D. Roekaerts'% F.J.J.Rosendal, F.G.J. Buitendijk, Ph. Daskopoulos, A.J.N. Vreenegoor and H. Wang Koninklijke/ShelI-Laboratorium, Amsterdam P.O. Box 38000, 1030 BN, Amsterdam, The Netherlands ~"~Also at Delft University of Technology, Faculty of Applied Physics, 2628 CJ Delft, The Netherlands Abstract Computational Fluid Dynamics (CFD) involves the numerical solution of conservation equations for mass, momentum and energy in a flow geometry of interest, together with additional sets of equations reflecting the problem at hand. In this paper the current capabilities of CFD for chemical reactor engineering are illustrated by considering a series of examples from industrial practice. These examples form the basis for the identification of future trends and potential stumbling blocks. Recent progress in the predictions of flow in baffled stirred tank reactors is reviewed. Improved turbulence modelling and high-performance computing have a particularly important part to play in the future applications of CFD to this type of reactor. In the next example, which concerns non-Newtonian, non-isothermal flow in an extruder, CFD turns out to be helpful in understanding the r61e of temperature profiles and residence time distribution in determining the product quality. In this application area limitations in available computing power and experimental validation of rheological models are both important. CFD simulations of gaseous turbulent flow with competing parallel and consecutive reactions constitute the third example. Apart from the need for accurate kinetic rates, the main problem in this case is the optimal choice of closure of chemical source terms in the mean transport equations for species mass fractions. Several models have been implemented in a commercial'CFD code and validated with plant data, thus allowing CFD model calculations to contribute to an improved reactor design. Finally, a discussion of future needs and expectations of the chemical processing industry with respect to CFD modelling, mainly focussing on multiphase flow, is given. 1. General Introduction Computational Fluid Dynamics (CFD) is a discipline that encompasses the numerical solution of the equations of motion (mass, momentum and energy) in a flow geometry of interest, together with subsidiary sets of equations reflecting the problem at hand. Three examples of such subsidiary equations are (1) equations for turbulence quantities in the context of the Reynolds-averaged formulation of the equations of fluid motion, (2) equations describing chemical species present in the flow and (3) equations describing the dynamics of solid particles, liquid droplets and gas bubbles dispersed in the flow. All these types of equations may well need to be considered in a reactor engineering problem. Further complicating factors could be that the flow occurs in a complex geometry, the fluid has non-Newtonian properties, and/or the chemical kinetics (homogeneous and heterogeneous) contains many steps. All these factors cannot always be taken into account in an appropriate way in a CFD calculation, which means that there is no universal answer to the question of what the rSle of Computational Fluid Dynamics should be in the development of new chemical reaction processes or the optimisation of existing ones. It is important to note that the term 'subsidiary' is used here to describe the transport, diffusion and reaction of chemical species only from a fluid flow perspective and is not, of course, intended to imply that they are of minor importance. In fact, they are of paramount importance, 1569 1570 C.K. HARRIS et al. as the purpose of chemical reactors in industry is the transformation of feedstocks into useful products. Indeed, from the point of view of a reaction engineering specialist it is the flow that is 'subsidiary' and, until relatively recently, it was always treated using a combination of highly idealised models. Although, in many cases, optimising the flow may only result in an increase in selectivity of one or two percent, the low margins in the industry mean that this can represent a staggering increase in profitability. The realisation of this has already motivated the application of detailed flow models to specific processes involving chemical reactions, while the continuing development of CFD technology, which has proceeded hand-in-hand with a corresponding development in computational capacity, has meant that CFD calculations are now increasingly being seen as an integral part of optimisation or development plans for a much wider range of chemical processes. These calculations help to identify the possible causes of problems occurring in existing systems and provide useful insights to be used in design of experiments. Furthermore, once they are accepted as a reasonably accurate description of the processes in a reactor they can be used for scale up. Here no attempt will be made to cover the past and projected achievements of CFD in all areas of chemical reactor engineering to which it has been applied or could be applied. We shall rather attempt to illustrate what the current capabilities of CFD modelling are by considering a series of examples from industrial practice. The examples will form the basis for the identification of future trends and potential stumbling blocks in the increasingly wide application of the CFD modelling approach, and its gradual emergence as a technology in its own right. In Section 2, we consider the application of CFD to baffled stirred tank reactors, and focus specifically on advances in treating the flow domain in these systems, which differs from virtually every other reactor type in that its shape changes with time. Simulations carried out using a variety of methods are presented, and compared with work carried out in the literature. Although an enormous amount of progress has been made in the past decade, simulation predictions for mean velocities and turbulence quantities are still not entirely satisfactory, even for single-phase non-reacting flows. Improved turbulence modelling and high-performance computing have a particularly important part to play in the future application of CFD to this type of reactor. Non-Newtonian, non-isothermal flow in an extruder forms the topic of Section 3. The main feature of this flow is the fact that the viscosity of the fluid, which is a novel type of polymer, exhibits a combined shear-rate and temperature dependence, as well as the fact that it is flowing through a complex geometry. In this application CFD modelling can help to optimise the process by improving the residence time distribution of the fluid so as to avoid hot spots which cause degradation of the polymer. The limits of current computing capability are again easily reached with this type of application, while another issue, particularly for fluids that are rheologically more complex than the one considered here, is experimental validation of the rheological models used. CFD simulations can provide the link between the rheological modelling and experimental results, thus enabling tuning of the models. The final example, described in Section 4, concerns gaseous turbulent flow in a reactor with competing parallel and consecutive reactions. Here the main physical problem to be considered is the influence of turbulent transport on the progress of the reactions. Various models for treating the turbulence-chemistry interaction are described, and their implementation into a commercially available flow solver is outlined. Again, currently available computer resources limit the application of the most sophisticated of these models to two dimensions. Comparison of the simulation predictions with experiments indicated that the widely used Eddy Breakup model yielded the worst predictions for the selectivities of the five product components considered. Computational fluid dynamics for chemical reactor engineering 1571 The present article is concluded, in Section 5, with a discussion of future needs and expectations of the chemical processing industry with respect to CFD modelling. These will mainly concern two-phase flow, as it is in this area that scale-up problems are most acute and general predictive tools are lacking. Space limitations preclude a detailed exposition, but the way we believe the development of CFD modelling should proceed in this fascinating and relatively uncharted area, and the pitfalls of some of the approaches that are currently appearing in the literature, will be outlined. 2. Computational Fluid Dynamics for Stirred Tank Reactors: New Methods and Old Problems 2.1 Introduction The application of Computational Fluid Dynamics to stirred tank reactors dates back to the late 1970's, but it is in the last few years that this particular reactor type has become a showcase for the development of CFD simulation technology. (3FD simulations of flow in baffled stirred tank reactors (BSTRs) are desirable because they can provide a useful supplement to the poorly established scale-up criteria that are traditionally used to design such reactors, in conjunction with the results of laboratory or pilot -scale tests. These criteria, as well as being of doubtful validity or even contradictory in some cases, involve global quantities such as power per unit mass, impeller Reynolds number and impeller tip speed, while local levels of the flow variables, particularly turbulence quantities which can vary by two or three orders of magnitude within the flow domain, may well have a crucial impact on reactor performance. The route from (3FD simulations of (single-phase) turbulent flow in BSTRs to the ultimate goal of predicting reactor performance usually requires phenomenological modelling approaches, backed up by extensive experimental validation, as a bridging step, though some multiphase (3FD simulations in BSTRs have been carried out recently. 2.2. Brief History of the Development of (3FD Technology as a Tool for Modellinq Turbulent Flow in BSTRs From the point of view of numerical simulation, a particular challenge presented by a baffled stirred tank is the fact that the shape of the flow domain (bounded by the liquid surface, the tank walls and baffles, and the impeller and its shaft) changes with time. This distinguishes it from most other reactor types in use in the chemical process industry. The first (3FD simulations of turbulent flow in BSTR's were two-dimensional, with the effect of the moving impeller blades and stationary baffles being smeared out in the tangential direction and represented as momentum sources and sinks, respectively, in the azimuthally averaged equation for the tangential mean velocity component. Appropriate momentum sources, depending on the type of impeller, were also included in the azimuthally averaged equations for the radial and axial components of mean velocity. Good examples of such simulations can be found in Refs. 1 & 2. The overall flow patterns obtained in these simulations were in encouragingly good agreement with experiment, and encouraged the further development of BSTR modelling using (3FD. Prior to the advent of (3FD modelling of BSTRs, the need for a more detailed description of the hydrodynamics than that implicit in the idealised (3STR model had long been realised in the industry. These included the network-of-zones model [3], based on experimentally observed flow patterns, while those without access to state-of-the art computing technology were attempting to work with grossly oversimplified versions of the hydrodynamic equations [4,5]. To the authors' knowledge, the first three-dimensional simulation of a BSTR to be reported in the open literature was the one by Middleton et a~ [6]. Prior to publishing their paper these authors, from 1(31and Imperial College, had already spent six years developing their code, which was specifically developed for BSTR modelling rather than for general-purpose 1572 C.K. HARRIS et al. applications. In addition to performing flow simulations, Middleton et al simulated a consecutive-competitive reaction sequence, with Iocalised injection of one of the components, so that the final outcome was dependent on the mixing produced by the flow circulation in the vessel• The simulations were able to provide a reasonably accurate prediction of the final yield, as a function of impeller speed, for a 30 L vessel and to demonstrate the failure of conventional scale-up criteria in extrapolating to a geometrically similar vessel with a capacity of 900 L. The reaction speed of the faster of the two reactions was, however, slow enough that only the large-scale features of the flow pattern could influence the progress of the reaction• The three-dimensional character of the simulation reported by Middleton et al did not extend to modelling the change in shape of the flow domain and, although the baffles were resolved, the effect of the impeller was smeared out in the azimuthal direction• It was to be several years more before CFD simulations of BSTRs were carried out that resolved both the baffles and the motion of the impeller blades. The intervening period saw the publication of increasingly large-scale simulations using 'impeller boundary conditions', in which the region swept out by the impeller is treated as a black box whose surfaces either contain inlets and outlets for the flow generated by the impeller, or which contains momentum sources distributed within it. Although some progress can be made by applying aerodynamic theory to the motion of an individual impeller blade [2], setting impeller boundary conditions normally requires experimental input, which is usually azimuthally averaged, and the resulting simulations smear out the effect of vortex shedding by the impeller blades and the subsequent interactions of the vortices with the baffles• Experiments have shown the latter effects to be quite striking [7]. These considerations motivated the development of alternative methods of treating the impeller which will be discussed at the end of the present section• 3D CFD simulations that have been carried out for BSTRs using impeller boundary conditions since the pioneering work of Middleton et al have involved both disc turbines [8-10] and downflow axial impellers [11,12] and have used in-house [8,11] as well as commercial [9,10,12] codes. Many of the simulations were carried out using 'standard' impeller and tank geometries. The work of Gosman et al [13] on Eulerian-Eulerian two-phase turbulent flow modelling in stirred vessels represents an important step forward, although the solids loading they considered was so low that there was negligible back-coupling to the fluid and the lack of coalescence and break-up mechanisms in their gas-liquid simulation work proved a serious shortcoming on attempting to compare the simulation work with experimental data. Returning to single-phase flow simulations, the comparison between 3D simulation predictions and experimental results for BSTRs reported in the literature exhibits generally good agreement for mean radial and axial components of mean velocity, while the tangential component is less well predicted, with the simulations often predicting spurious flow reversal near the tank floor and liquid surface. It is this aspect of the simulation predictions that is most sensitive to the choice of turbulence model, with anisotropic turbulence models yielding superior results• The only attempts to take the flow-induced deformation of the liquid surface into account have been for laminar flow - see for example Ref 9. Although some of the work reported in the literature yielded good results for the prediction of turbulence quantities in the impeller stream, elsewhere in the vessel predictions of turbulence quantities are generally poor • This is an area of major concern, as a knowledge of turbulence quantities, especially the turbulent dissipation ~, is needed in order to attempt to quantify processes on which the operation of commercial BSTRs operating in the turbulent flow regime depend, namely, the degree of micromixing, bubble break-up and turbulent diffusion. Reliable predictions of ~, in particular, would be particularly valuable as experimental values are always inferred rather than measured, since direct measurement would involve the simulataneous determination of nine independent velocity gradients at each measuring point (see, Eg., Ref. 14). Recent attempts at simulating consecutive-competitive reaction sequences, where turbulent Computational fluid dynamics for chemical reactor engineering 1573 micromixing limits the conversion rate of the fastest of the reactions, are reported in Refs. 15 & 16. After fitting coefficients in the turbulent micromixing models used to the data, moderate agreement with experimental work was obtained. The disadvantages of setting effective, steady, impeller boundary conditions, namely the neglect of transient coherent structures such as eddies shed by the trailing edges of the impeller blades as they move through the fluid, and the need for experimental input have already been mentioned. Compilations of the latter exist for standard impeller types, but a substantial experimental effort would be involved in providing impeller boundary conditions for, say, testing a series of novel impeller designs. An elegant way around this problem that does not involve any additional computational effort is to carry out CFD simulations in a frame that is rotating with the impeller, so that the latter can be explicitly modelled, while incorporating a momentum source model for the baffles in the rotating frame. This approach was applied by Hutchings etal. [17] for a BSTR and by Dong e t a l [ 1 8 ] for an unbaffled tank, in which case the boundary condition at the tank wall is simply set by imposing a tangential wall velocity in the rotating frame. For a BSTR, it is a natural step to azimuthally average the flow solution, and use the resulting values on the surface of the impeller swept volume to generate synthetic impeller boundary conditions, which can then be used as input for a conventional simulation in a frame of reference fixed with respect to the tank wall. The inner-outer method takes the approach just outlined one step further, by defining a fictitious cylindrical boundary with a radius intermediate between that of impeller blade tips and the baffle edges furthest away from the wall, azimuthally averaging the flow solution over this surface, and using the resulting values as an outer boundary condition for a second simulation in the rotating frame. Repeating this iteration process a number of times will result (hopefully) in a converged steady-state flow solution. This technique was applied successfully by Brucato et al [19]: the inner and outer boundaries used are depicted schematically in Figure 2.1a. Alternatively a single, common boundary could be used for both the inner (rotating reference frame) and outer (fixed reference frame) simulations, as illustrated in Figure 2.1b. A largescale flow simulation of this type, using a parallel version of STAR-CD, was reported very Common boundary recently by Jones [20]. for inner/outer method or sliding mesh simulations Outer boundary for inner zone simulation baffle baffle I Inner boundary for outer zone /~ simulation r (a) . . . . . I IL . . . . . I I I Impeller swept volume I (b) Figure 2. l a&b. Alternative subdivisions of flow domain into zones for the inner-outer method. Finally, in the sliding mesh method, one again has the arrangement of meshes depicted in Figure 2.1b, but the simulation is carried out in time-dependent fashion with the relative motion of the meshes and consequences for data transfer across their common boundary taken into account. There are a number of specialised techniques available for doing this. Luo et al [21] report a simulation of this type, though one can infer from their article that the flow is laminar rather than turbulent. 1574 C.K. HARRISet al. In the following section we report in-house simulations using three different commercial c o d e s and all three methods of treating the flow domain. The simulation results are compared with the experimental data of Ranade & Joshi for a Rushton turbine [22] and a pitched-blade impeller [23]. We consider this data the most comprehensive available in the open literature. 2.3. Review of In-House Simulation Work 2.3.1 Simulations Using Impeller Boundary Conditions These simulations were carried out using FLUENT for the tank fitted with a Rushton turbine and CFDS-FLOW3D for the tank fitted with a pitched-blade impeller. The data reported by Ranade & Joshi in Refs. 22 & 23 was used to set the impeller boundary conditions. Both cases involved standard geometries, with the height of the liquid zone being equal to the tank diameter, and the tip-to-tip diameter of the impeller and baffle width being one third and one tenth, respectively, of the tank diameter. The tanks were fitted with four baffles, equally spaced around the tank wall, and the impellers were of a standard type. The Rushton turbine was located in the centre of the liquid zone, while the pitched-blade impeller was placed one third of the way from the tank floor to the liquid surface, along the axis of the vessel. Turning first to the simulations of the tank fitted with the Rushton turbine, the liquid surface was modelled as a no-slip boundary in order to be able to consider the impeller centre plane as a symmetry plane. This would be appropriate for a tank that was completely filled with liquid and operated with a lid: it is not entirely clear whether this was the case in the experiments performed by Ranade & Joshi. In any case, provided that one limits comparisons with experiments to the bottom two-thirds of the tank, the effect of the boundary condition at the liquid surface will be fairly minor. By utilising the symmetry plane just mentioned, together with the four-fold rotational symmetry exhibited by the flow domain with azimuthally-averaged impeller boundary conditions, it is in principle sufficient to discretise one-eighth of the domain in order to perform flow simulations. Because of difficulties in implementing 90-degree periodicity in FLUENT, however, a 180-degree sector of the lower half of the tank was simulated. The grid size used was 20 × 20 x 23 in the axial x radial x tangential directions, amounting to a 40 x 20 x 46 grid for the entire tank. These dimensions compare with an effective grid size for the entire flow domain of 60 × 15 x 32 as used by Ranade & Joshi in Ref. 8. Simulations were carried out using the standard k-~ turbulence model, an RNG version implemented in FLUENT, and the standard differential Reynolds stress model. Figure 2.2a shows radial profiles of mean axial velocity (measured vertically downwards) in a mid-baffle plane at various distances (0.933, 0.833, 0.667 and 0.5 in units of tank radii) below the impeller centre plane, while Figure 2.2b depicts corresponding flow profiles for the turbulent kinetic energy nearer to the tank centre (0.2, 0.333 and 0.5 radii below the impeller centre plane). The k-£ and RNG turbulence models yield good agreement with the experimental data for mean axial velocity, except very close to the tank floor, where they underpredict it. The Reynolds stress model, on the other hand, severely underpredicts the mean axial velocity throughout the bottom quarter of the tank. The standard k-~ model overpredicts the turbulent kinetic energy in the wall jet below the impeller, with the RNG and Reynolds stress models yielding better results. However, the Reynolds stress model underpredicts k in the bottom quarter of the tank (not shown, though this trend is already seen for the bottom curve in Figure 2.2b). Figure 2.3 shows radial profiles of tangential velocity a little way below the impeller and near the tank floor. In contrast to the experimental data, the k-~ and RNG turbulence models predict a fluid rotation in the opposite sense to the impeller rotation throughout much of the flow domain. This deficiency of the simulations is largely corrected by using the Reynolds stress model, but quantitative agreement with the experimental results is generally poor. Computational fluid dynamics for chemical reactor engineering O r'----"--"l ~3( 1575 IZlR=0'2 I 0 0 0 0 0 0 o o ~ i!~ x ~°lo. o o o o o o o o,b,~'~ ,°1g ~0~ 4434 k-'i-'- (a) ....... DIMENSIONLESS RADIAL COORDINATE, r / R Figure 2.2. I. • 1.0 ooO°! (b) DIMENSIONLESSRADIAL COORDINATE, r / R Radial profiles of (a) mean axial velocity and (b) turbulent kinetic energy in a baffled stirred vessel fitted with a disc turbine. Grey solid line, black stippled line and grey stippled line: in-house simulations performed using the standard k-~, RNG k-~ and Reynolds stress turbulence models, respectively. (a) Bullets: Measurements Ref. 22, Solid line: k-~ simulations Ref. 8 - viceversa in (b). Dashed line in (a) and open circles in (b) are k-~ simulations from Ref. 8 performed with the vessel wall treated as a free surface. I z,R= o.2'¸':'i................... - .......................... ..................................... !iiil,t o~, ~"-°"I I UJ i z/R=0.933 Z I Z t • 0 , ~ ~ m . ' , ~ ~- N +0"ZO Figure 2.3. 0-8 0.t 04 0,Z 0'~ 0"4 0",~ 0.6 DIMENSIONLESS RADIAL COORDINATE, r / R 0.9 1.5 Profiles of mean tangential velocity in a baffled stirred vessel fitted with a disc turbine. Notation is as in Figure 2.2a. 1576 C.K. HARRISet al. In order to test the effect of grid refinement, especially on the predictions of turbulence quantities, additional simulations were carried out, using the standard k-£ turbulence model, on a 38 x 38 x 44 grid (corresponding to a 76 x 38 x 88 discretisation of the entire flow domain). Some of the results for mean velocities and turbulence quantities are compared with coarsegrid simulation results in Figures 2.8 and 2.9. The main conclusions to be drawn from these results is that, while predictions of mean velocity profiles exhibit little change, overprediction of k in the wall jet worsens on refining the grid. Nearer to the tank axis, and especially in the bottom quarter of the tank (not shown), the refined grid yields improved predictions of k. Some simulations were also carried out with the liquid surface modelled as a free surface instead of a no-slip boundary. These confirmed, as expected, that the different boundary condition had little effect on the predictions of mean velocities and turbulence quantities in the bottom half of the tank. The simulations of the tank fitted with a pitched-blade impeller were carried out on a 90-degree sector using CFDS-FLOW3D. The standard k-~ and differential Reynolds stress models were used. The grid size used was 60 x 58 x 16, which amounts to 60 x 58 x 64 for the full tank. Ranade and Joshi used a grid size corresponding to the entire flow domain of 45 x 30 x 16 [11]. Figure 2.4a shows radial profiles of mean radial velocity in the mid-baffle plane at a sequence of increasing distances measured downwards from the impeller centre plane. The top profile in Figure 2.4a is located just above the impeller swept volume, while the remaining three are at increasing distances below it. Figure 2.4b shows corresponding profiles for turbulent kinetic energy. Except, curiously, for the profile that is the smallest distance downstream of the bottom face of the impeller, where impeller boundary conditions were set, the simulation predictions for the mean radial velocities are in good agreement with the experimental data. The two turbulence models yield very similar predictions in this region of the tank. The striking feature of Figure 2.4b is the remarkable underprediction, by an order of magnitude, of the turbulent kinetic energy in the impeller stream. This is not reflected in Ranade and Joshi's simulation results, obtained using their own code, and merits further investigation. The predicted mean tangential velocity is, unlike the case of the tank fitted with the disc turbine, in fair to good agreement with experiment below the impeller, with the two turbulence models yielding very similar results. This situation changes drastically on approaching the liquid surface. Figures 2.5a and 2.5b compare horizontal sections just below the liquid surface through the velocity field predicted by the k-£ and Reynolds stress models respectively. The k-~ model predicts that most of the fluid is swirling in the opposite direction to the impeller, while the Reynolds stress model predicts a seemingly more realistic flow pattern. Unfortunately Ranade & Joshi did not report any measurements in the upper part of the tank. 2.3.2 Simulations Performed using the Inner-Outer Method For the tank fitted with the Rushton turbine, these simulations were carried out using the code STAR-CD, as part of an evaluation of the code, and used the standard k-£ model. A common inner/outer boundary was used, which was located midway between the tips of the impeller blades and the inner edges of the baffles. Although, in principle, symmetry would allow one eighth of the outer zone to be simulated and one twelfth of the inner zone, for ease of programming the simulations were in fact carried out using a quarter of each zone, corresponding to the lower combined symmetry of the inner and outer zone. The grid size covering the entire flow domain was 66 x 31 x 50. The sensitivity of the results to a local refinement of the impeller region was also tested (STAR-CD allows embedded meshing). Figure 2.6a compares simulation results for mean axial velocity at the same locations as Figure 2.2a and should be compared with the latter. Axial velocities are underpredicted, suggesting the smearing out of the wall jet by numerical diffusion. The turbulence in the wall jet is grossly underpredicted, which tends to corroborate this view, and Figure 2.6b shows that serious underprediction of the level of turbulence even extends to the impeller stream (these 1577 Computational fluid dynamics for chemical reactor engineering % 7:> x (..) ,st" o ,< z m Z _o 09 ° == ~ o~ t' i z / R = 0.356 ~.................... I.-tel ._1 Z o_ ~-'~ ~--~... . .* . . ~. . o |oo ° o ~ ~ o - 3 o ~ ~ 09 ooeo 3~ (a~ ,, ! DIMENSIONLESS RADIAL COORDINATE, r / R (b) J 0.I , O'Z ~ 0"3 DIMENSIONLESS i ° I 0.~ , i O'Z z/R=0.533 I i e..e ~3-; ! ,:,,~ RADIAL COORDINATE, 0,9 r / R Figure 2.4. Radial profiles of (a) mean radial velocity and (b) turbulent kinetic energy for a baffled stirred vessel fitted with a pitched blade impeller. Grey lines: in-house simulations performed with the standard k-E model using impeller boundary conditions (solid line) and the inner-outer method (stippled line). Stippled black line: in-house simulations performed using the differential Reynolds stress model (impeller boundary conditions). Solid lines: (a) k-~ simulation predictions from Ref. 11 and (b) experimental measurements reported in Ref. 23. Symbols: (a) experimental data reported in Ref. 23 and (b) simulation predictions from Ref. 11 (bullets: k-~ simulation predictions with wall functions, and open circles: k-E simulation predictions with tank wall considered as a free surface). \ \ \ (a) Figure 2.5 (b) Horizontal cross-section through simulated flow field in baffled stirred vessel fitted with pitched blade impeller, 1 mm below liquid surface. (a) k-~ turbulence model, (b) differential Reynolds stress turbulence model. C. K. HARRIS et al. 1578 o x i I i >.' I- > x / r-R~ DtMENSIONLESS RADIAL COORDINATE <-- (a) z o O3 Figure 2.6. CI C.-) 0,2 0.3 0.~, O.S ~,s 0-7 o:e g"9 Radial profiles of (a) axial velocity and (b) turbulent kinetic energy in a baffled stirred vessel fitted with a disc turbine. The profiles shown in (b) are located in the impeller centre plane. Stippled grey line: Present simulation results using the innerouter method (k-~). Bullets (a) measurements reported in Ref. 22 and (b) k-E simulation predictions reported in Ref. 8. Solid line and dash-dotted line: measurements reported in Ref. 22. Dashed line and open circles represent k-~ simulations carried out in Ref. 8 with the tank wall treated as a free surface. i,o DIMENSIONLESS RADIAL COORDINATE, r / R Axial Velocity Radial Velocity Tangential Velocity Turb. Kinetic Energy Turb. Dissipation/100 • . . . - ...... + ...... ........ B ................ x Axial Velocity ~ ~ ......Y ....... Radial Velocity . . . . . .......... Tangential Velocity ...... .......... " Turb. Kinetic Energy ........... ................... Turb. Dissipation/100 ........ . . . . . 2 ~> i i " [] [] [] []..G...~..~.-m--6"' 1 • .~,.o.-- ~:~.::.:.:.:.:.:L .iLL ....................................... ~'g.:g.g...~'"~'"~'"'~'"';~ x x x x x x x × .~ Z . . 6 . L I - L ~ - ~ . . * - . A _ . , ~ - . ~ . ~ _ I L ~ - . ~ . ~ _ & - ~ g ~ ~ ............. -1 (a) Figure 2.7. (b) -2 0.01 0.015 0.02 0.025 0.03 0.035 radius 0.04 0.045 0.05 0.055 (a) Surface of inner mesh used in inner-outer simulations of a baffled stirred vessel fitted with a pitched blade impeller. (b) Comparison of impeller boundary conditions on the bottom face of the impeller swept volume (lines) with corresponding values recovered after iterating between the meshes. Computational fluid dynamics for chemical reactor engineering results can 2.9 below). predictions depicted in 1579 be compared with the ones obtained using impeller boundary conditions in Figure Again, the direction of swirl of the fluid below the impeller is wrongly predicted: the obtained using the present method are broadly similar to the k-~ simulation results Figure 2.3. CFDS-FLOW3D was used to perform simulations of the tank fitted with the pitched-blade turbine. Overlapping inner and outer zones were employed in the simulations, as depicted in Figure 2.1a. The standard k-~ model was used. An impression of the impeller geometry, as modelled in the CFDS-FLOW3D code, can be obtained by inspecting Figure 2.7a, which shows surface meshes on the impeller blades, hub and shaft and on the outer boundary of the inner zone. These simulations were limited to a test of the inner/outer method in which a converged flow solution obtained using impeller boundary conditions was used to set the boundary conditions for the inner zone simulation, and the resulting flow field was azimuthally averaged over the surface of the impeller swept volume. The flow variables on the bottom face of the impeller swept volume were then compared with the impeller boundary conditions applied in the original simulation. The results of this comparison are shown in Figure 2.7b. While the mean velocity profiles are recovered with a moderate degree of accuracy, the profile of turbulent kinetic energy is underpredicted by a factor of two (though it has the correct shape) and the turbulent dissipation is grossly underpredicted, being a factor of about 30 too small. It is interesting to note that the inner zone simulation yields values of k in much closer agreement with experiment below the impeller than the simulation using impeller boundary conditions (see Figure 2.4b). 2.3.3 Simulations Performed using the Sliding Mesh Method The recently released FLUENT sliding mesh module was used to perform simulations, using the k-~. model, of the tank fitted with the Rushton turbine. Figures 2.8a and 2.8b, respectively, compare simulation results with experimental data for radial profiles of mean axial velocity and turbulent kinetic energy below the impeller. The locations of the profiles are as in Figures 2.2. The profiles of mean axial velocity predicted using the sliding mesh method are similar to those obtained using impeller boundary conditions, while the turbulent kinetic energy profiles exhibit severe underprediction. This underprediction is, however, not as drastic as that obtained using the inner-outer method (not shown). Figure 2.9 shows that the severity of the underprediction of k in the impeller stream is comparable to that obtained using the inner-outer method (cf. Figure 2.6b). The simulation predictions obtained using impeller boundary conditions and the standard k-£ model on coarse and fine grids are also shown in Figure 2.9 for comparison. Finally, the prediction of mean tangential velocities exhibits the same deficiencies as for the k-£ model implemented using impeller boundary conditions. 2.4 Discussion The fact that the CFD simulations of flow patterns in BSTRs reviewed in the preceding section were carried out using off-the-shelf commercial codes is of particular relevance to industry, since most industries involved in CFD modelling have chosen to purchase such codes rather than develop in-house software. A serious disadvantage of these codes is that the details of the algorithms used to obtain a flow solution are often only partly available to the user. Thus, on considering the simulation results obtained with impeller boundary conditions, it is difficult to determine whether a difference in the quality of the prediction of turbulent kinetic energy in the impeller stream obtained using different codes (compare Figures 2.2b and 2.4b) is due to a different way of treating the inlet boundary condition. For the simulations performed using the inner-outer and sliding grid methods, it is probably safe to conclude that insufficient spatial resolution, especially near the impeller, is responsible for the low levels of turbulence observed in all cases. A striking feature of the study reported in the previous section was the excessive computational run time involved - as long as 20-50 1580 C.K. HARRIS et al. o 0 0 0 0 0 >_I-- $° o > ~0 u) Z~ ~ Z 2 0 0 00 0 0 0 ~_ lO o_ 0 GO m "~ -.&,~._...J . . . . . . . ,.&__ I , I @# # ~ . , i ' ~0~0 O ' C l~ C 0 O D ~ C; C~ +" , , , ........ L DIMENSIONLESS RADIAL COORDINATE, r / R Figure 2.8. G (b) DIMENSIONLESS RADIAL COORDINATE, r / R Radial profiles of (a) mean axial velocity and (b) turbulent kinetic energy in a baffled stirred vessel fitted with a disc turbine. Thick grey and black solid lines: in-house simulations performed using the standard k-~ model and impeller boundary conditions on a coarse and fine grid respectively. Grey stippled line: present simulation results carried out using the sliding mesh method on a coarse grid (k-~). Remaining lines and symbols pertain to measurements and simulations carried out in Refs. 8 & 22 and are as in Figure 2.2. e / R ~ C,.O~ x ! ~,.1 w o z w 0,0b~ N t 0 ~ ~.0,L w G 0[ .-J.--~.~--L--- D I M E N S I O N L E S S RADIAL C O O R D I N A T E Figure 2.9. 0-1L 0,g r-R. Radial profiles of turbulent kinetic energy in the impeller centre plane for a baffled stirred vessel fitted with a disc turbine. Symbols as in Figure 2.8b, with the dash-dotted line corresponding to the thin solid line in the latter. Computational fluid dynamics for chemical reactor engineering 1581 CPU hours, on an RS6000 580 workstation, for the FLUENT sliding mesh simulation, with the inner-outer method requiring a factor of four less computing time. Parallel versions of all the most commonly used commercial codes have recently or shortly will become available and the next year or two will see these being applied to perform impeller-resolved BSTR simulations on ever larger grids. The work of Jones et al [20], which involved an effective grid size for the entire fluid zone of almost two million cells, provides a good foretaste of future trends in this area. Another very recent development is the application of Large Eddy Simulation using the lattice Boltzmann method to this flow [24]. Very fine meshes of up to 14 million cells have been used and the passage of the impeller blades has been modelled by applying appropriate timedependent constraints on the fluid velocities within the impeller swept volume. Good results for flow quantities, including turbulent intensities, have been obtained. The traditional method of using impeller boundary conditions will still be used for routine applications, especially if these are extended to include multiphase, reacting flows, for some time to come. Here the emphasis will be on turbulence modelling, especially on the development of an appropriate multiphase flow model that is able to take account of relevant features of the phase dynamics in a physically sensible way. Even for single-phase hydrodynamics, certain features of the flow patterns are wrongly predicted that can be attributed to deficiencies in the turbulence model used. For example, poor predictions of the magnitude and even the sign of the mean tangential velocity are a recurring feature of the simulation results reported in the previous section. It is for the reasons mentioned in this paragraph and the previous one, in conjunction with the central importance of the BSTR within the chemical process industry, that the stirred tank reactor will remain an exacting benchmark carried along at the forefront of CFD simulation technology. 3. Numerical Simulation of Polymer Flow through an 8-to-O section 3.1 Introduction The design of commercial scale extruders is often based on experiences built up in the past. Although these experiences are very important, they are difficult to use in, for example, the scaling up of existing extruders to larger machines. Experiences are also often related to a certain polymer which makes it difficult to use them for the design of a set-up for a new, totally different product. To obtain more insight in the physical phenomena occurring in an extruder detailed measurements might be very useful. However, the results of the measurements can not always fully explain the observed phenomena. In that case, numerical simulations can give additional information. The present example concerns a measurement programme carried out on a commercial scale double screw extruder. Although a lot of information was obtained from the measurement results, some questions remained. Numerical flow simulations on this complex geometry could be performed at KSLA because of the availability of an SP2-machine which contains 30 RS6000/390 nodes. The calculations were carried out on one "fat" node with an internal memory of 1024 Mbytes and a swap space of 1500 Mbytes. These extensive numerical simulations resulted in plausible explanations for the observed phenomena, thus proving again the strong combination of experiments together with numerical simulations. The experiments performed consisted of measurements of temperatures on a commercial scale double screw extruder to study the distribution of the polymer melt temperature over the die plate. In order to obtain the temperature data the die plate has been fitted with 24 thermocouples distributed evenly over the die plate area. The results of these measurements are depicted in Figure 3.1. The polymer melt temperature distribution has been determined by interpolating and extrapolating the 24 measured values. Figure 3.1a shows the temperature of the polymer melt at the start of the extrusion where a nearly homogeneous temperature can be seen of about 245 °C with a variation of only a few degrees Celsius. Figure 3.1b shows the temperature distribution near the end of the extrusion. The homogeneous temperature 1582 C.K. HARRISet al. distribution at the start has disappeared and two "hot spots" can be distinguished at 4 and 10 o' clock, with temperatures up to 265 °C, while the other areas measure 255 - 260 °C. During the measurements degraded material flowed occasionally from specific die holes at 4 and 10 o'clock positions (seen from the downstream side of the die plate). It was unlikely that the degradation of this material was caused by a temperature of only 5 % higher. In order to try to identify the reason for the outflow of degraded material at the two "hot spots" the measurements were supplemented with numerical simulations performed using the commercial POLYFLOW code. 3.2 Setup of the CFD problem The simulations involved non-isothermal flow of a generalised Newtonian fluid through the 8to-O section of a co-rotating twin screw extruder. The quantities computed were streamlines, and temperature and residence time distributions. The 8-to-O section of a double screw extruder is the location where the two screws end and the 8-shaped barrel gradually turns into a cylinder. The cylindrical part guides the flow towards the die plate. The simulations have been performed on a geometry which included some previously proposed modifications to the original set-up. As a result, the geometry used for the numerical simulations differed slightly from the geometry used during the experiments. This may affect the detailed comparison of the present simulations with the experimental observations. The computational mesh and flow domain are shown in Figure 3.2a. The mesh includes the tips of two co-rotating screws. The present state-of-the-art software does not allow the user to account for a rotating flight on both screws. The rotation of the tips induces the highest shear rates on the rotating tip walls. The mesh has been refined towards the tip walls as a consequence. In order to calculate the viscosity q(y ,T) a power law was used to describe the shear-thinning behaviour of the polymer. The POLYFLOW code has been modified and extended to account for the non-standard temperature behaviour of the polymer. The boundaries are displayed in Figure 3.2a. The boundary where the polymer melt is in contact with the outer steel wall is associated with zero slip and zero heat flux boundary conditions. Outflow conditions hold on the exit of the cylindrical section where developed flow has been assumed as a consequence. Developed flow is also assumed at inflow together with a homogeneous inlet temperature. The screw speed was set at the same value as during the measurements. A zero heat flux has also been assumed along the walls of the rotating tips. 3.3 Numerical aspects The five equations of motion expressing the conservation of mass, momentum and energy must be solved on the computational mesh for the appropriate material parameters and boundary conditions. Viscous heating will be taken into account while inertia and gravity are neglected. Starting from a trivial solution POLYFLOW will not converge as a result of the highly non-linear nature of the problem. Therefore, a special numerical strategy must be followed, to be defined by the user. Boundary conditions and material parameters generate the following orders of magnitude for the dimensionless numbers of interest: Re = O(10-4), Pe = O(104), Ga = O(106). The order of magnitude for Reynolds number Re indicates that, as expected, viscous forces dominate over inertia forces. The P6clet number Pe shows that convective heat fluxes dominate the diffusive fluxes and the Galilei number Ga demonstrates that gravitational forces may be neglected with respect to viscous forces. The high Peclet number emphasises the convective (hyperbolic) nature of the energy equation, giving rise to non-physical oscillations in the temperature field caused by numerical instabilities if this equation is not treated in a special Computational fluid dynamics for chemical reactor engineering (a) Figure 3.1. 1583 Co) Temperature distribution of polymer melt at die plate (a) just after exhaustion commences and (b) 408 seconds later. flow domain / l no slip q=0 [] outflow ,Z2- , ~; nflow \, f ~ '%.,~e.~.. -~ -,.:. :. -. rotating screw tips inflow streamlines I :~i ¸ I IC) • ,~ ~ i :J:~.,,,: n ~ |_,,,,. -2~o.o Figure 3.2. 100 x heat conduction Simulation results showing (a) Computation geometry, (b) Flow streamlines of fluid subjected to maximum viscous heating (i.e. originating between the screws), (c) Temperature distribution over flow domain boundary and (d) Distribution of fluid whose residence time within the flow domain is greater than 25 seconds. Computational fluid dynamics for chemical reactor engineering 1585 way. The following strategy has partially proven to be successful in achieving a well converged solution free of numerical wiggles: - - - set the heat conduction coefficient k at 100 times its actual value (decrease Pe), compute an isothermal solution by means of a double evolution scheme on a gradually increasing flow rate and a gradually increasing screw speed, with this solution, start a non-isothermal calculation with an evolution scheme that gradually decreases k by a factor of 100 to its original value. The results of the non-isothermal simulations associated with the evolution on the heat conduction coefficient k showed that convergence worsened with decreasing k (increasing Pe). At higher Pe number, the convective term in the energy equation starts to dominate, and numerical oscillations appeared in the solution of the temperature field. In six steps k could be reduced by a factor of 4.5 (a factor of 100 was required). In the corresponding converged solution numerical oscillations were clearly present. The oscillations are generated on the rotating screw tip walls at the location where the tips change from cylindrical into conical. Mesh refinement around this location produced some improvement but could not eliminate the problem. Note that the present mesh results in a job requiring 1300 Mbytes RAM, challenging the limits of "fat node" the computer system used; the mesh cannot be refined much further. Because of the numerical oscillations present in the solution after six evolution steps, we have decided to present the results of the first non-isothermal solution, associated with a heat conduction coefficient that is a factor 100 too high. 3.4 Results of Isothermal Simulations The 3D flow simulation has been visualised with the Data Visualizer software. Although stationary flow is considered, a clear insight into how the computed quantities vary as the flow proceeds through the geometry is greatly assisted by animation techniques, involving various rotations of the 3D domain and the release of massless particles to visualise streamlines. A set of animations has been composed into a video-presentation of the flow. Some snapshots of the animations are discussed below. We recall that the present results are obtained with the appropriate material data, except for the heat conduction coefficient k which is multiplied by a factor 100. Figure 3.2b shows some streamlines. Observed from the downstream side both screws perform a clockwise rotation. When massless particles are released at inflow from the centre of the location of maximum viscous heating (between the screws) they eventually flow from both tips to arrive at a 2 and 8 o'clock position in the cylindrical section downstream. Figure 3.2c shows the temperature distribution across the boundaries of the 8-to-O section. The temperature increase is small (1 °C) as a result of the high heat conduction coefficient; decreasing k will generate higher temperatures confined to smaller volumes. Nevertheless, it is clearly visible how hot material on the screw tips is convected to the 2 and 8 o'clock positions mentioned above. An inhomogeneous temperature distribution across the cylindrical portion of the geometry downstream of the tips occurs as a result. Note that the present 8-to-O section has a different cylinder diameter than the 8-to-O section of the experimental set-up. Evidently, this will affect the convection of hot polymer from the screw tips to the die plate which can explain the difference of an angle of 60 ° for the calculated and measured position of the two hot-spots. In addition to the correct prediction of two "hot-spots" occurring at the die plate, numerical calculation of the residence time showed another interesting result: The volume of fluid having a residence time that exceeds 25 s is displayed in Figure 3.2d. The regions of high temperature and high residence time obviously coincide which is unfavourable for a polymer. These numerical results can be related to observations during the measurements where degraded material occasionally flowed from die holes at 4 and 10 o'clock positions on the downstream side of the die plate. Moreover, temperature measurements have shown hot spots at these locations of the die plate. 1586 C.K. HARRISet al. Initially, the simulations were planned to calculate the convection of hot polymer from the screw tips downstream to the cone and flow channel towards the die plate with the purpose to investigate the location of hot spots in a cross-section where temperatures actually had been measured in practice. The mesh was extended accordingly. A flow simulation could not be performed, however, as a result of the excessive RAM requirements generated by this mesh (note that the 840-0 section alone already requires 1300 Mbytes RAM). This simulation may become feasible with the installation of the parallel beta-version of POLYFLOW on the SP2. 3.5 Conclusions 3D numerical simulations have given a clear insight into the flow of a polymer melt through the 8-to-O section of a co-rotating double screw extruder. The results predict an inhomogeneous temperature distribution at the outflow that has also been measured in practice. The coinciding regions of high temperature and high residence time could explain the occasional outflow of degraded polymer at two specific positions. Only qualitative results have been obtained for the temperature field as a result of the occurrence of numerical oscillations at low values for the heat conduction coefficient. The nature of the oscillations is presently under investigation. 4. Turbulence-Chemistry Interaction in a Gas Phase System 4.1 Introduction We consider a single phase system where reactants here called A and B are mixed in a nozzle and subsequently react exothermally to form C and D. Apart from the main reaction leading to the desired product there are other parallel and consecutive reactions leading to unwanted byproducts E, F, G and H. The eight species participate in five reactions: A+B---~C+D A+B--+F A+B-+C+G D+B---~C+E A+D-->C+H The molecules in this system are built from three different elements. The species have mass fractions Y=(Y1 ..... Y8) and the source term for species o~(in kg m -3 s-1) is pSr~with p the density. Both p and S depend on __Y,the temperature Tand the reference pressure. To one species in the model correspond several species in reality, but one model species always represents molecules with a unique element composition. The rates of the five reactions have been determined in laboratory scale experiments and have been fitted with expressions which take into account the nature of the underlying detailed kinetic mechanism. In previous work a simpler version of the system was considered [25,26]. There continues to be considerable incentive to improve the selectivity for this system and a sequence of activities has been defined involving model calculations, cold flow tests and verification in a plant, leading to an improved reactor design. The objective of the project is an improvement in the weight % of D in the mixture of D, E, F, G and H by say 5 %. Starting from the observation that the model for turbulence-chemistry interaction is an important component of the CFD model for the reactor we have made a comparative study of the performance of three turbulence-chemistry interaction models: the Mean Value (MV) model in which mean reaction rates are expressed in terms of mean concentrations and mean temperature, the extended Eddy Breakup (EBU+) model in which the MV reaction rate is replaced by a turbulence driven reaction rate when chemical processes are faster than turbulent mixing processes and the probability density function (PDF) model in which the exact Computational fluid dynamics for chemical reactor engineering 1587 definition of mean reaction rates in terms of the joint PDF of all independent thermochemical variables is used and an approximation of the PDF is sought for. In this work the shape of the joint PDF of the chemical variables is assumed, with the first and second moments to be determined from modelled transport equations [27,28]. The more general Monte Carlo solution technique in which the Lagrangian evolution of notional fluid particles is computed to represent the PDF [25,26,29,30] has not been used here because general purpose computer programs for the latter approach are not available yet for the case of complex three dimensional geometries• The assumed shape PDF model we used is a simplified one in which the fluctuations in composition and temperature due to turbulence are represented by a joint PDF for two variables. The first variable is a conserved scalar which is a measure of the local equivalence ratio of the mixture (the fraction of the mass coming from one of the two inlet streams).The second variable is a reaction progress variable which is a measure of the conversion reached locally in the reactor. Validation is done by comparing the predicted yields with the measured yields for existing reactors. 4.2 Models for Turbulence-Chemistry Interaction The transport equations for velocity, mass fractions and enthalpy are solved numerically. As usual for turbulent flow at high Reynolds number, the quantities solved for are not the instantaneous values, but rather the mean values, averaged over timescales much greater than that pertaining to the most persistent turbulent fluctuations. Because an exact closed system of equations for mean values does not exist, modelled equations have to be used. Here we use the standard k-~ model for the closure of the turbulent velocity field. The transport equation for the mean mass fraction of species ct is As usual for variable density systems, density weighted averages are used in the modelling [29]. The gradient diffusion assumption is used for the closure of the turbulent transport of thermo-chemical variables. To close the mean chemical reaction rate three models have been used: • Mean value (MV) closure consists in the assumption that the mean reaction rate is obtained by substituting the mean mass fractions and the mean temperature in the kinetic expressions. It is a good model in the limit that the reactions are slow compared to the turbulent mixing rate and the fluctuations are small. • The extended Eddy Break-up (EBU+) model is the same as the MV model except for the fact that when any reaction rate tends to be very fast the kinetic expression is replaced by a phenomenological rate proportional to the frequency associated with the large turbulent eddies. For example when SD.1 is the kinetic rate for formation of D by the first reaction: ILl ~ l),1 ] with n.J - CF.SU W ~ • min , % 1588 C.K. HARRISet aL where W~ is the molar mass of species o~.The eddy breakup conversion is proportional to the turbulent mixing frequency, i.e. turbulent dissipation rate ~ divided by turbulent energy k, and to the concentration of the limiting reactant. CESUis an empirical coefficient. • Probability density function (PDF) closure uses the fact that the mean source terms can be written exactly as Here _~ and 74 are possible (instantaneous) values of the mass fractions and the temperature and f is their joint probability density function, which is to be determined. Closing the equations involves the use of an approximation for the PDF. The main challenge in the assumed shape PDF method is to make reasonable assumptions concerning the joint PDF. In the most general case all species and enthalpy would be fluctuating independently and many moment equations would have to be solved. In practice the fluctuations in these variables are correlated. First of all there are a set of rigid constraints following from element conservation: because of conservation of elements in chemical reactions the element mass fractions Zk, k=1,2,3, which are linear combinations of the mass fractions satisfy a transport equation with vanishing chemical source term. In the presence of only two inlet streams and assuming equal diffusivity for all species, a unique normalised nonreacting scalar can be defined as: Z k -Zk, 1 z~-z~,~ where the suffices 1 and 2 denote the values in the first and second inlets respectively. If we assume that the syslem is adiabatic, the total specific enthalpy is a linear function of ~, and hence also a conserved scalar. These observations reduce the number of independent fluctuations from eight to five. However, this is still more complex than necessary. Because of the correlations following from the structure of the reaction system (e.g. all reactions rates respond similarly to a change in temperature) further simplification is possible. The simple model we propose, following related literature on combustion problems [27], is that the essential characteristics of the fluctuations in the space of independent scalar variables can be captured by considering just two fluctuating variables, one linked to the mixing of the two inlet streams (the conserved scalar ), and one related to the conversion (a progress variable ). The reaction progress variable is zero before reaction starts and one in the final state. In general several alternative definitions are possible. We choose: Y. y;;~tx (~, Yih,7,,.,,d,tct, ) Here YDmax is the maximal attainable mass fraction of D, given a certain value of the conserved scalar and given certain mass fractions of the unwanted byproducts (E,F,G and H). In this two variable approach, the PDF closure consists in the assumption that the mean source terms are given by I ] 0 () Computational fluid dynamics for chemical reactor engineering 1589 Here _Y~and -/4 are the instantaneous value of the mass fractions and the temperature corresponding to a possible value ~t of the mixture fraction and a possible value r~of the progress variable and .f~,ris the joint PDF of ~ and r. To characterise the PDF in the usual way one would have to solve the equations for the variance and covariance of ~ and r in addition to the equations for their mean values. Now and rare assumed to be statistically independent so that their covariance is identically zero. In fact, the main reason to introduce the scaled variable ris that ~ and rare more likely to be independent than ~ and Yo. To obtain instantaneous values of the mass fractions needed in the evaluation of the source term some extra hypothesis is needed. Our hypothesis is that fluctuations in mass fractions of the byproducts can be neglected, i.e. we assume that for species E,F,G and H t h e instantaneous mass fractions are equal to the mean mass fractions. It then follows that fluctuations in progress variable are strictly related to the conversion of the first reaction, producing D. Because the fluctuations have to respect energy conservation, fluctuations in progress variable are also directly connected to fluctuations in temperature. Attributing fluctuations in progress of the reactions completely to the first reaction, is not the same as assuming a mean value closure for all other reactions. There can be a change in all mean reaction rates through the difference in predicted mean temperature.For given values of conserved scalar and progress variable the mass fractions of the species occuring in the first reaction A,B,C and D can be obtained using element conservation and solving a linear system of equations.. We have Ygt=f I YDmax (~1.,y~ byproducts ), where YDmax is an important auxiliary variable. Once ~/ is known, 74 can be obtained from the caloric equation of state. Using element conservation and statistical independence of ~ and r it can be shown that it is sufficient to solve for the mean mass fractions, the variance of a the conserved scalar and the variance of the mass fraction of species D, to determine the mean and variance of ~ and r a n d by assump~ion, the PDF. 4.3 Results The commercial CFD code CFDS-FLOW3D, extended with user defined subroutines incorporating the three proposed models, has been used to calculate the turbulent reacting flow in a cylindrical jet stirred reactor similar to the one described in Refs. 25,26 for which plant data were available at a base case operating condition and at conditions with increased pressure and with increased throughput and changes in preheat temperature. In the computations with the PDF-model a very simple assumed shape for the PDF, the double delta PDF, was used both for mixture fraction and for progress variable. The results of the model calculations have been compared to plant data and the following observations have been made: The total selectivity to D is predicted better by the MV model and the PDF model than the EBU+ model. In particular the EBU+ model gives incorrect trends as the pressure increases. It predicts a 2.5 % decrease in D in favour of all byproducts. This is not seen in the plant data which show only 1 % decrease predominantly in favour of F. In fact, the EBU reaction rate SAEBu is the same for the first three reactions. Because the kinetic rates of these reactions are different and the EBU rates are equal the point at which these reactions become mixing limited is different. This is in agreement with the ideas underlying the model, but it is counterintuitive that parallel reactions between the same two species are treated differently with respect to mixing limitations. It is also clear that the precise point of switching between kinetic rate and EBU rate, which is a phenomenological property expressed in the model constant, is interfering directly with the prediction of selectivity. The reaction rate of the EBU model scales 1590 C.K. HARRISet al. linearly with pressure since turbulence frequency is rather independent of pressure and the rate is of first order in the reactant concentrations. On the other hand, to leading order, the kinetic rates scale quadratically with pressure. Consequently at higher pressure the number of reactions for which the EBU reaction rate is used and also the spatial domain where it is used becomes larger and the extent to which the kinetic rates are taken into account is correspondingly reduced. This explains why the model is sensitive to changes in pressure and fails to predict the measured trend. The differences between the MV model and the PDF model are very small. It is in agreement with the expected behaviour of the kinetic rates that the PDF model predicts a higher fraction of D than the MV model, but the predicted difference is only about 0.2 %. The differences in selectivity to D between the different operating conditions considered are of the order 1 % and are predicted equally well by both models. The extra information provided by the PDF model has been studied in the form of contour plots of the relevant variables (Mean and variance of mixture fraction, mean and variance of mass fraction of species and of progress variable, standard deviation of temperature). It is seen that fluctuations in mixture fraction and in progress variable occur in different regions of space. Significant fluctuations in mixture fraction occur inside and close to the nozzle only (Region I in Figure 4.1 ). A perfectly mixing nozzle would lead to a uniform mean mixture fraction with value following from the incoming fluxes of A and B. This perfectly mixed state is not reached inside the nozzle and a certain unmixedness remains in region I. Significant fluctuations in progress variable occur in the center of the reaction at the border of the region where the D-formation rate is high (Region II in Figure 4.1). A+BI in[ L Figure 4. 1. Schematic diagram of reactor showing zones of unmixedness (I) and reaction (11). Both fluctuations in mixture fraction and in progress variable lead to fluctuations in temperature. In the former case this is because the two inlet streams are at a different temperature, in the latter case this is because states with different status of progress have a different temperature. Considering the strength of the fluctuations and the residence time in the different parts of the reactor it can be understood why the fluctuations have only a small impact on the predicted selectivity. 4.4 Conclusions The observations of the model validations have lead to the recommendation that for the subsequent modelling studies of alternative reactor geometries the predictions of the EBU+ model should not be used in the comparison and ranking of different reactor geometries. The MV model and PDF model can both be used in comparative studies of promising reactor geometries, but since the computational cost of the MV model is smaller and the difference in predictions for the reactors studied until now was found to be very small it was recommended to use the MV model with the refinement of the PDF model to be used as a check at the end. Computational fluid dynamics for chemical reactor engineering 159 I 5. Discussion and Possible Future Trends in CFD for Chemical Reactor Engineering Two themes common to the three applications of CFD treated in the preceding sections are the fact that they were all carried out using commercial codes developed and marketed by third parties and they were all limited in some way by the computer resources currently available. In the case of the application to stirred tank reactors, this limitation meant that simulations using the inner-outer and sliding-mesh methods had to be restricted to relatively coarse grids, with dire consequences for the prediction of turbulence quantities. For the extruder application, computer memory was the limiting factor, and meant that it was not possible to refine the modelling by extending the grid. Finally, in the application to turbulent reacting flow, expected computer run-time limitations precluded the implementation of the most sophisticated PDF modelling approaches currently available. All developers of commercial CFD codes either already have or are in the process of developing parallel versions of their software that are suitable for running on workstation clusters with both shared and non-shared memory. Many large industrial organisations have such clusters and the future will see a large increase in the use of parallel CFD to solve industrial problems. The vast majority of CFD modelling for practical problems in chemical engineering (meaning complex geometries and/or complex physics) will continue to be carried out using commercial codes. This is because it is generally not considered cost-effective for industrial users to develop CFD simulation codes themselves, even though such an effort could be tailor-made to the users' applications rather than, as with an 'off-the-peg' general-purpose code, being liable to lack some detailed features required to treat some specialised applications while containing other unnecessary ones. It is our experience that the documentation provided by commercial code developers is not detailed enough to answer many specific issues which arise naturally in the process of setting up a CFD simulation, especially those of a numerical nature. Solving numerical problems takes up a considerable amount of time, especially when attempting largescale CFD modelling in complex geometries. Commercial code developers are continually introducing new physical models into their software but the effect of implementing the models on the numerical performance of the code is an issue that deserves more attention. Finally, all three applications treated in this paper involved single-phase flow simulations. In many types of reactor, catalyst particles and/or gas bubbles are present in the fluid. Two examples that can be singled out as being of particular importance to our industry are riser reactors, in which the major issue is the design of internals and inlet devices to alleviate the maldistribution of the catalyst particles across the riser cross section, so increasing selectivity, and slurry reactors, whose commercial implementation has to overcome notorious scale-up problems. In a CFD simulation the simplest way of representing the dynamics of these particles or bubbles is to treat them as structureless objects experiencing a variety of forces, of which the most important are gravity/buoyancy and hydrodynamic drag. For low enough disperse-phase hold-ups, the influence of the particles or bubbles on the continuous fluid phase and on each other may be neglected. In this case, the velocity field and turbulence properties of the continuous phase may first be obtained by solving the equations of fluid motion for that phase in the absence of particles or bubbles, and then solving for the latter individually as a postprocess, introducing them into the flow domain with spatial and velocity distributions appropriate to the process concerned. However many processes of interest, including the two examples given above, involve dense dispersions of bubbles and/or particles. The interactions of the bubbles and particles on the fluid can be taken into account, at vastly increased computational expense, by solving the fluid flow equations in the presence of a distribution of forces representing the hydrodynamic interactions between the disperse phase and the continuous phase. Various ways can be envisaged to treat particle-particle interactions, the 1592 C.K. HARRISet al. simplest probably being a modification of the hydrodynamic drag term to include a dependence on the solids hold-up. In the case of bubble-bubble interactions, break-up and coalescence are additional processes that must be treated. Furthermore, for the large bubbles that arise in fiuidised beds and bubble columns, the evidence is mounting that the detailed dynamics of bubble shape plays a crucial r61e in the hydrodynamics of the entire system. As an alternative to solving subsidiary equations in which particles or bubbles are tracked through the continuous-phase flow field, the proposal has been made to treat these phases as if they were continuous and to model them on an equal footing with the continuous phase [31]. This approach, termed Eulerian-Eulerian multiphase flew modelling, is considered an attractive alternative, for dense dispersed multiphase flows, to the Eulerian-Lagrangian approach outlined in the previous paragraph because it is computationally much more efficient: however, there remain serious problems concerning the physical interpretation of a disperse phase modelled as a continuum analogue as well as the appropriate relations to assign to effective properties of the disperse phase such as the viscosity. Recently there has been a flurry of articles in the literature applying Eulerian-Euterian CFD modelling to bubble and slurry columns: unfortunately, in those cases for which quantitative comparison with experimental results is made, the authors either use empirical input for the radial distribution of the gas hold-up, or treat a coefficient multiplying the lift force as a free adjustable parameter. It is also unfortunate that some authors are rather circumspect about which parameters have been fitted to experimental data (compare, for example, Refs. 32 and 33). The best way forward may well be a judicious combination of Lagrangian and Eulerian techniques such as described in Ref. 34. If an embedded modelling approach is adopted, validation of the various submodels employed is essential if the end result is to capture the complex interplay between the various physical phenomena involved in the dynamics of a multiphase reactor in a meaningful and useful way. References 1, 2. 3. 4. 5. 6. . . 9. 10. Harvey, P.S. & Greaves, M., (1982), turbulent flow in an agitated vessel. Parts I & II. Trans. Inst. Chem. Eng. 60, 195-200 & 201-210. Pericleous, K.A. & Patel, M.K., (1987), The source-sink approach in the modelling of stirred reactors. Phys. Chem. Hydrodynam. 9,279. Mann, R., (1986), Gas-liquid stirred vessel mixers: Towards a unified theory based on networks-of-zones. Chem. Eng. Res. & Des. 64, 23-34. Fort, I., Obeid, A. & Brezina, V., (1982), Coll. Czech Chem. Commun. 47,226. Fort, I., (1983), Turbulent flow in an agitated vessel. Comment & reply. Chem. Eng. Res. Des. 61,136-137. Middleton, J.C., Pierce, F. & Lynch, P.M., (1986), Computations of flow fields and complex reaction yield in turbulent stirred reactors, and comparison with experimental data. Chem. Eng. Res. & Des. 64, 18-22. Van der Molen, K. & Van Maanen, H.R.E., (1978), Laser Doppler measurements of the turbulent flow in stirred vessels to establish scaling rules. Chem. Eng. Sci. 33, 11611168. Ranade, V.V. & Joshi, J.B., (1990), Flow generated by a disc turbine: Part II. Mathematical modelling and comparison with experimental data. Trans. IChemE 68A 34-50. Brucato, A., Ciofalo, M., Grisafi, F. & Rizzuti, L., (1990), Computer simulation of turbulent flow in baffled and unbaffled tanks stirred by radial impellers. In 'Computer Applications to Batch Processes', Cengio, Italy, held 28th-30th March, 1990. Kresta, S.M. & Wood, P.E., (1991 ), Prediction of three-dimensional turbulent flow in stirred tanks. AIChE J. 37, 448-460. Computational fluid dynamics for chemical reactor engineering 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 1593 Ranade, V.V., Ooshi, J.B. & Marante, A.G., (1989), Flow generated by pitched blade turbines I1: Simulation using k-~ model. Chem. Eng. Commun. 81,225-248. Bakker, A., (1992), Hydrodynamics of stirred gas-liquid dispersions. PhD Thesis, Delft University of Technology, The Netherlands. Ch. 3. Gosman, A.D., Lekakou, C., Politis, S., Issa, R.I. & Looney, M.K., (1992), Multidimensional modelling of turbulent two-phase flow in stirred vessels. AIChE J. 38, 1946-1956. Kresta, S.M. & Wood, P.E., (1993), The flow field produced by a pitched blade turbine: Characterisation of the turbulence and estimation of the dissipation rate. Chem. Eng. Sci. 48, 1761-1774. Ranade, V.V. and Bourne, J.R., (1991), Reactive mixing in an agitated tank. Chem. Eng. Commun. 99, 33. Bakker, A. & Fasano, J.B., (1993), Time dependent, turbulent mixing and chemical reaction in stirred tanks. Preprint, Paper No. 70C presented at the AIChE Annual Meeting, St. Louis, Mo., November 1993. Hutchings, B.J., Weetman, R.J. & Patel, B.R., (1989), Computation of flow fields in mixing tanks with experimental verification. Paper TN-481, presented at the ASME Meeting held San Francisco, 1989. Dong, L., Johansen, T., & Engh, T.A., (1994), Flow induced by an impeller in an unbaffled tank- I1. Numerical modelling. Chem. Eng. Sci. 49,3511-3518. Brucato, A., Ciofalo, M., Grisafi, F. & Micale, G., (1994), Complete numerical simulation of flow fields in baffled stirred vessels: The inner-outer approach. Proceedings, IChemE Symposium Series No. 136, 155-162. Jones, D., Dimirdzic, I., Krishna, R. & Robinson, D., (1995), Use of parallel CFD for demanding chemical process applications. In 'Chemputers Europe I1', Proceedings of the conference and exhibition held Noordwijk, The Netherlands, 10th-12th October, 1995. Published by Chemical Engineering. Luo, J.V., Gosman, A.D., Issa, R.I., Middleton, J.C. & Fitzgerald, M.K., (1993), Full flow field computation of mixing in baffled stirred reactors. In Proceedings of the t993 IChemE Research Event, held Birmingham, U.K., 6th-7th January, 1993. Published by the Institution of Chemical Engineers, Rugby, U.K., pp. 657-659. Ranade, V.V. & Joshi, J.B., (1990), Flow generated by a disc turbine: Part I. Experimental. Trans. IChemE 68A, 19-33. Ranade, V.V. & Joshi, J.B., (1989), Flow generated by pitched blade turbine I: Measurements using laser Doppler anemometer. Chem. Eng. Commun. 81,197-224. Eggels, J.G.M., (1995), Direct and large-eddy simulation of turbulent fluid flow using the lattice-Boltzmann scheme. Jubilee issue of J. Heat & Fluid Flow, to be published. Roekaerts, D., (1991), Use of a Monte Carlo PDF method in a study of the influence of turbulent fluctuations on selectivity in a jet-stirred reactor. Applied Scientific Research 48,271-300. Roekaerts, D., (1992), Monte Carlo PDF method for turbulent reacting flow in a jetstirred reactor. Computers and Fluids 21, 97-108 Correa, S.M. & Shyy, W., (1987), Computational models for continuous gaseous turbulent combustion, Progr. Energy Combust. Sci. 13, 249-292. Bockhorn, H., (1988), Finite chemical reaction rate and local equilibrium effects in turbulent hydrogen-air diffusion flames. Twenty-Second Symposium (International) on Combustion, the Combustion Institute, 655-664. Pope, S.B., (1985), PDF methods for turbulent reactive flows. Prog. Energy Combust. Sci. 11, 119-192. Pipino, M. & Fox, R.O., (1994), Reactive mixing in a tubular jet reactor: a comparison of PDF simulations with experimental data. Chem. Eng. Sci. 24B, 5229-5241. Ishii, M., (1975), Thermo-fluid Dynamic Theory of Two-Phase Flow. Eyrolles, Paris. 1594 32. 33. 34. C. K. HARRISet aL Torvik, R. & Swendsen, H.F., (1990), Modelling of slurry reactors: a fundamental approach. Chem. Eng. Sci. 45, 2325-2332. Wachi, S. & Yates, J.G., (1991), Comments on modelling of slurry reactors - a fundamental approach. Comment and reply by Torvik & Swendsen. Chem. Eng. Sci. 46, 1528-1530. Lapin, A. & Lebbert, A., (1994), Numerical simulation of the dynamics of two-phase gas-liquid flows in bubble columns. Chem. Eng. Sci. 49, 3661-3674. Acknowledgements - The authors wish to thank Dr. Gert Colenbrander for a critical reading of the manuscript and the copyright holders of Refs. 8 (Gordon and Breach Science Publishers SA) and Refs. 11 & 23 (The Institution of Chemical Engineers) for their kind permission to reproduce figures from these articles in the present work. List of Symbols A,B, .... H D R s~ T Utip Y~, z, f k k q r r W Z symbols of chemical species laminar diffusivity radius of stirred tank reactor chemical source term of species cc temperature component of velocity vector impeller tip velocity mass fraction of species c~ mass fraction of element k probability density function turbulent kinetic energy (not Section 3) heat conduction coefficient (Section 3 only) heat flux radial coordinate (Section 2) progress variable (Section 4) spatial coordinate mean axial velocity vertical coordinate in stirred tank reactor Greek letters y ~1 p (~ shear rate turbulent dissipation dynamic viscosity density mixture fraction turbulent Schmidt number Symbols denoting averages in turbulent flow Reynolds average density weighted average