CFD Based Qualification of Mixing Efficiency of Stirred Vessels

In this work, we focus on the most crucial units in a chemical technology, the chemical reactors. Using a commercially available CFD software package, COMSOL Multiphysics, 3D mathematical models of a batch reactor with different impeller geometries have been investigated. The reasonable agreement between the experimental and simulation results indicates the validity of the developed CFD model. The effect of the impeller design, e. g. number of blades on the mixing efficiency is evaluated based on the simulation studies. The proposed measure to determine the energy efficiency of mixing (i. e. mixing index) is based on the calculated velocity field and energy usage. The information about the homogeneity of the mixed phase in the system can be extracted from the developed velocity field. Hence, we proposed histograms of velocity fluctuations on a logarithmic scale as an efficient tool to measure the achieved homogeneity of the phase in case of different impellers and rotational speeds.


Introduction
During the last few decades, computational fluid dynamics (CFD) method has become a very powerful tool in the process industry primarily to analyse the hydrodynamic behavior. In general, CFD methods rely on the numerical solution of the Navier-Stokes equations, thus flow velocity, velocity magnitudes, and pressure values can be calculated in the geometric space investigated. A particularly significant benefit of the CFD simulators is that they can include entire apparatuses in three dimensions, regarding either single phase or multiphase systems [1,2]. This property of the CFD simulators makes them useful in a large number of research areas, including the study of mixing phenomena. Fluid mixing is a fundamental operation in the chemical process industry and the stirred tank or autoclave is a commonly used reactor type in many kinds of chemical technologies; therefore these systems are well described. Regardless, when it comes to design work, simple experimental correlations are still commonly used, e.g. Reynolds number, Froude number, power number or the Euler number which is derived from the previous ones. This approach has its limitations when it comes to the design of stirred vessels with complex internal structure.
On the other hand, in today's mixed equipment multiple impellers and/or baffles are built in to reduce the dead zones and intensify the rate of processes take place in the reactor. CFD can be used to model multi impeller systems, e.g. gas dispersers and biochemical reactors [3]. The effect of impeller and baffle configurations on gas dispersion can be evaluated and the optimal structure of the disperser can be determined. CFD can be applied to simulate more complex systems, such as fluidized bed reactors in detail e.g. fast pyrolysis of biomass in fluidized bed reactors [4]. It can also be used to study process intensification problems [5].
Fluid mixing is usually characterized by tracking the dissipation of an indicator component injected into the bulk of the fluid [6]. One of the earliest attempts to introduce a mixing number belongs to Danckwerts [7]. He proposed two largely independent processes to describe fluid mixing. The breaking-up of the fluid reduces the size of the clumps, while the interdiffusion tends to equalize the concentration differences between the neighbouring regions of the mixture. Tumasz and Thiffeault created a "topological braiding exponent" to characterize the complexity of the braid generated by the motion of three or more fluid particles [8]. Another method for quantifying mixing is MixNorm, which has been created by Mathew et al. [9]. This method roots in the concept of weak convergence that averages over all scales of mixing to produce a single number characterizing the extent to which an advected scalar field has mixed. In addition, several methods have been developed to quantify mixing (Poincaré maps, periodic point analysis, Lyapunov exponents, etc.) [10].
Fluid mixing takes place by a combination of mechanisms: convection caused by bulk motion (stirring or macromixing), stretching, folding, and diffusion [11]. In the stirred liquid there can be certain areas called dead spots where only diffusion takes place. In general, the more the kinetic energy of the liquid is, the less is the share of the dead spots from the fluid volume, therefore the mixing time will reduce. Therefore, mixing time can be directly calculated from the velocity of the mixed fluid [12]. The velocity can be measured by various indirect techniques such as the Laser Doppler Anemometry (LDA) [13,14], Particle Image Velocimetry (PIV) [15,16], or the less common Positron Emission Particle Tracking (PEPT) [17,18]. Furthermore, it can also be calculated with a proper CFD simulator.
Despite the fact that mixing efficiency of stirred vessels can be characterized by the evolving velocity field, there is no common method which can be easily and widely used in comparison various velocity fields gained by the application of different types of impellers or different rotation speeds. In this paper, we present a method based on the distribution, mean, and standard deviation of the velocity magnitude, a method which can be easily applied in the comparison of various velocity profiles from CFD simulations.
Another important field to be addressed for qualification of mixing efficiency is the power consumption that affects the operating costs directly; therefore, it is an important variable in designing of stirred tanks. Power consumption is often determined by the power number of the agitator. The dimensionless number based power usage calculation has its limitations; theoretically, it can only be applied for simple impeller geometries (e.g. Rushton turbines [20]). In case of a more complex agitator, investigating the effect of the dimensions of the applied agitator on the power number has its limits because of the high number of potential variables affecting the power consumption. Detailed and time-consuming experimental investigations are required to find the parameters for the recommended empirical equations.
Power number can be determined without such difficulties from the power consumption. Iranshahi et al. determined the power number of the agitator experimentally in a specific system using Maxblend impeller and pointed out that the value depends on the number of baffles, therefore, the geometry of the vessel [21]. However, finding the exact correlation between the power number and the geometry requires detailed inspection on that matter. CFD simulation can be an excellent and cost-effective tool to calculate mixing efficiency using CFD simulators even with complex impeller geometries [22][23][24].
Despite the large amount of work on this topic, there is no consensus on a single measure of mixing that can be used for different processes nor is there a single mixing index that can be easily implemented for different scenarios from numerical simulations to laboratory experiments to industrial processes. In this work, mixing power requirement of different pitched blade and Rushtonturbines were calculated from the force awakened on the impeller surface at different revolution speeds and temperatures. We also propose a new method to correlate the calculated power requirement to the measured values.

Experimental and modelling methods 2.1 Experimental system
A small-scale laboratory vessel (cylindrical design with 100 mm diameter and 1 dm 3 volume) was used to investigate the power consumption of different impellers. The agitator has a variable speed drive (IKA Eurostar Power-Control Visc mixing motor). Different impellers were manufactured by a 3D prototype printer based on CAD drawings. The overall layout of the system can be seen in Fig. 1 (a). To measure the power intake of the drive, we used a Voltcraft Energy Check 3000 type electric consumption meter. Measurement of power intake of Rushton turbines with different numbers of blades (shown in Fig. 1 (b)-(e)) and four-bladed blade turbine (shown in Fig. 1 (f)) was investigated at different rotational speeds (from 50 min -1 up to 400 min -1 ). The power requirement of the system was measured at different revolution speeds and different type of impellers. A mixture of ethylene-glycol and propylene-glycol (7:3) was used in our study as its above mentioned physical properties -viscosity and density -show high temperature-dependency (shown in Fig. 2). Note that while the experimental system is intrinsically homogenous, the proposed method of velocity field fluctuation evaluation can be applied to any one-phased system, therefore apparent mixing level of the liquid phase in the vessel can be concluded from our results as well.
The distance between the motor and the reactor were as small as possible to ensure that the vibration in the system is minimal. Every experiment was repeated multiple times to ensure the reliability of our measurements.

CFD model of the stirred vessel
The Rotating Machinery model of COMSOL Multiphysics was used to calculate the flow pattern of the vessel, in which the fluid motion caused by impeller rotation is described using the Navier-Stokes equations and the RNS k-ε turbulence model [25]. The model is robust and widely used in many CFD applications. However, solving the model requires a finer grid near walls of the equipment, where no-slip condition (i.e. velocity is 0) is defined.
To model the rotation of the shaft and the impeller, the built-in Moving Mesh interface of COMSOL Multiphysics was used which defines the shape-change of the geometry (mesh) as the displacement of the spatial frame relative to the material frame [26]. The impeller and the shaft boundaries were treated as rotating walls and the rotational speed was determined as a function of the space coordinates. COMSOL Multiphysics built-in solver was used (PARDISO), developed by Schenk et al. [27].
A three dimensional geometric model of the stirred vessel was built because the impellers can only be described using spatial coordinates. This geometry was then divided into a large number of finite volume elements (mesh). The quantity of the mesh elements was determined in order to reach the required computational accuracy in a manageable simulation time. The results of the mesh independence study are shown in Appendix 1. In Fig. 3. (a) there is an example of the used mesh in case of a three-bladed Rushton turbine, which contains 63036 tetrahedral mesh-elements. The elements were refined near the edges of the rotating walls. This applies in general to all models built. Further increasing the number of mesh elements does not have a significant effect on the results, hence we considered that as optimal and applied it in . We used the same order of magnitude mesh elements to solve each model. The developed time-dependent model was solved for 10 s because at the end of this time interval a relatively stationary (quasi-stationary) velocity field is evolved in the experimental apparatus. To reach the solution easier a step function was used to achieve the specified rotational speed in each simulation. In this way, the rotational speed is not jumped up from zero to the predefined value, but it takes 1 s to reach it. In order to extract the required information from the calculated velocity, a 10 x 10 x 10 grid was described within the geometry in the program as axes were divided into equidistant parts as shown in Fig. 3. (b). In each point of the grid, the local velocity value was collected. The size of the grid was chosen as so to give a fair representation of the developed velocity magnitudes while also keeping the memory usage of the calculations with this three-dimensional matrix manageable. Fig. 3 (c) shows the developed velocity field inside the reactor equipped with a four-bladed Rushton turbine with 400 min -1 revolution speed and Fig. 3 (d) shows the velocity magnitudes in grid points, showing 10 s is enough to achieve a quasi-stationary state. The impeller rotational speed varied between 50 and 400 min -1 with 50 min -1 increments both in the simulational and experimental environment. The effect of the variable density and viscosity on the power intake of the drive was investigated between 20 °C and 90 °C with 10 °C steps.
Velocity magnitudes extracted in the grid points were plotted in histograms in case of every impeller and rotational speed. The logarithmic abscissa of these histograms represents the fluctuation of velocity values, and the ordinate represents the frequency of these fluctuations in the grid. We investigated the power requirement of the impellers at every rotational speed and temperature. Power number and mixing time for these impeller geometries were introduced in our previous work [28].

Results and discussion
To obtain a reliable picture of the homogeneity of the emerged velocity field (i.e. the deviation of velocity) in the mixing apparatus, it is not enough to examine the last moment of the process. Therefore, we calculated how much time the impellers require to turn around twice, and the calculated velocity fluctuations during this time interval were plotted on histograms. In order to easier interpret of the histograms, mean values (MV) and normal standard deviations (NSTD) were calculated for different impellers at the investigated temperature and rotational speed combinations.
where n refers to the number of the grid points were used which is a vector in our case because the defined grid points rearranged as a vector after the simulation. If the NSTD value shifts to higher values, that means higher shear rate between fluid elements and results in better mixing; this can be derived from the width of the histograms. The effect of the rotational speed is shown in Fig. 4-6. At 20 °C (Fig. 4), it can be seen that at higher rotational speed the mean values rise at every impeller types, so the fluid can be exchanged quickly between individual fluid elements. The higher the average fluctuation and the smaller the standard deviation are, the higher the level of mixing is in the apparatus. If the rotational speed is increased eightfold, the value of MV increases two orders of magnitude.
On the basis of the results in Fig. 4, it can be stated that the rotational speed has a significant impact on the temporal fluctuation of velocities in the device. However, there is no trend in the values of NSTD, so it can be stated that most of the velocity fluctuation in fluid elements are slightly different from each other. Fig. 5 shows the results of the homogeneity of the velocity field in the mixing apparatus at 50 °C. Because the viscosity of the mixed fluid phases decreases with increasing the temperature, velocity fluctuations can develop in a lower range. Hence, at this temperature higher mixing level can be achieved than at lower temperatures at the investigated rotational speeds. At 50 °C, the NSTD values are extremely high using five-and six-bladed Rushton turbines at 400 min -1 . In this case, the fluctuation is very large, so is the mixing level resulted in case of these turbines at low revolution speed. At 400 min -1 the four-bladed blade turbine results similar mixing level to three-and four-bladed Rushton turbines. Furthermore, the effect of the rotational speed on the MV values decreases, and with the increasing of the rotational speed at all impeller types the average of the velocity fluctuations increases. Fig. 6 shows the results at the highest investigated temperature (90 °C), i.e. at the lowest viscosity. The trend, which has been experienced earlier, remains in case of each impeller type, so the average fluctuation of the velocity values increases with the growth of the rotational speed, but the NSTD values of these fluctuations decrease. The rotational speed has the greatest impact on the fluctuations at this temperature that is with the growth of the rotational speed the fluctuation can increase by four orders of degree. The electric power requirement for an impeller depends largely on the shape and the size of the impeller and of course on the rotational speed. To define a mixing efficiency measure we should know how much energy should be used to reach a specific level of mixing. Hence, we calculated the power requirement for each impeller type. From Eq. (4) the consumption can be determined from simulation data. The force awakened on the surface of the blades is calculated with a surface integral, which is parametrized as a function of x and y coordinates in COMSOL Multiphysics. By multiplying this force by the length of the lever arm, in this case, the length of the mixing blades (r), a torque type result is obtained which then is multiplied by the angular velocity (ω) to get the power requirement of the impeller at the specific conditions.
The electrical energy requirement is measured at 20 °C for each impeller in the rotational speed range of 50-400 min -1 with the help of the experimental system. For each impeller type we performed three parallel measurements, and we used the average of these. We measured the power requirement of the shaft (with no impeller attached to it), and the combined power requirement of the shaft and the impeller. The difference of these values gives the measured electrical power requirement of the impeller. To compare this with the results of simulations, we divided the measured data with the simulation data. This gives a ratio type value (i mix ) which provides the link between the measured and the calculated values. The reason why the measured and calculated energy requirements are different is that only the applied mixing motor power requirement can be measured and we do not have information about the characteristic of the motor. For a better understanding of this ratio Eq. (5) was defined, in which the properties of the mixed fluid are present, therefore they can be connected to our factor. Parameters of Eq. (5) were determined with a global search algorithm, particle swarm optimization in MATLAB program [29]. With this algorithm, we minimalized the relative error (square of the difference of i mix calculated by Eq. (5) and divided by the i mix calculated from measured and calculated data) and we fitted three parameters to calculate the ratio. With this method, we can calculate the electrical power requirement based on simulation results at higher temperatures.
In Eq. (5) Re represents the dimensionless mixing Reynolds-number, while m is the number of blades. Table 1 shows the identified parameter values. Fig. 9 shows the measured and calculated electrical power requirements in the dependence of the mixing Reynolds-number at 20 °C. It can be seen that the measured and calculated power requirement values give a good approximation for most impeller types. There is a major difference in the case of the four-bladed Rushton turbine. This can be originated from some kind of systematic error during the measurement, which will be checked in further investigations.    proposed mixing efficiency measures the higher value is the better since it means that high mixing level can be achieved in the mixing apparatus with low energy input. It can be seen that the energy-specific fluctuation of velocity values decreases with the increase of the rotational speed at each temperature, in case of each impeller type. In other words, for achieving better mixing the required power will increase exponentially which corresponds to the classical relationship. From these results, the degree of homogeneity for the given mixing problem can be found out.
The power required to obtain the necessary mixing level can be determined too, and it can be stated out that whether it is sufficient or not. In general, we can say that the energy-specific velocity fluctuations in the case of four-bladed blade turbine are located between three-and four-bladed Rushton turbines.

Conclusions
Computational Fluid Dynamics is widely used to study the mixing performance of different types of impellers or to design new or improved equipment. The fact that the application of CFD modelling can greatly reduce the amount of experimental work required makes its application highly advantageous. In our work, mixing efficiency of various impellers was studied quantitatively. Different operational parameters, namely the effect of rotational speed (50-400 min -1 ), number of blades and material properties of the mixed phase (density, viscosity in the temperature range 20-90 °C) were investigated in case of three-, four-, five-, six-bladed Rushton turbines and a four-bladed blade turbine. Velocity magnitudes extracted in grid points were applied to calculate velocity fluctuations plotted in histograms in case of every impeller and rotational speed to get a picture on homogeneity. The electric power requirement of the impeller was determined and it was validated by measurements on the real system. The described methods are highly flexible, it can be suggested that any impeller geometry can be studied in any one-phase fluid; therefore it is a method of high performance for impeller selection, equipment design and process intensification in industrial applications as well. Future work might include the investigation of mixing efficiency in multiphase flows, extending the applicability of these methods to a number of other areas with industrial relevance such as fermentation processes.