Numerical Analysis of Rotational Speed Impact on Mixing Process in a Horizontal Twin-Shaft Paddle Batch Mixer with Non-Newtonian Fluid

Nowadays, the priority of food industry processes is based on mixing of the various flavour pastes with the air. Optimal mixing process in the results means the instantaneous preparation of the homogenous and easily spreadable paste while consuming the lowest power. For achieving the brand precise and repeatable flavour, of major importance are, beside the filling parameters, the mixing setup parameters. This work presents the influence of rotation velocity, as one of the setup data, on: the paste quality and the mixing index, obtained by the CFD. The paste was prepared of soy flour and water, and described by the experimentally determined rheological characteristics. Mixing vessel, of the twin-shaft paddle type, introduced mixing by the laminar shear-thinning. This mixing process was modelled through 3D fully transient numerical simulations employing the volume of fluid (VOF) multiphase model of the Eulerian–Eulerian approach. In two cases, the influence of the angular velocity to the mixing-index value can be neglected. At the mixing area two regions were marked off. First one and dominant, was the region of active mixing preoccupied with high shear flow, described with high mixing index. Opposite, in the passive mixing region, the flow was almost rotational, represented by low mixing indexes. The rotation velocity change made the significant influence to wall shear stress levels achieving the highest values at the tips of the paddles.


INTRODUCTION
Dispersive mixing of highly viscous, non-Newtonian fluids are one of the most common and important operations in the food industry.When such materials are involved, the mixing process is generally in the laminar region, and full attention must be paid to the design of the mixing elements in order to obtain the efficient mixer with desired performances.Engineers in the chemical industry recognized the opportunities of the numerical modelling and simulation in obtaining reliable numerical results that lead to efficient processes and design solutions of complex requirements [1][2][3][4][5][6], simultaneously reducing the number of validation tests and costs.Food science and engineering traditionally have been the basis for the development of the process technology related to the manufacture of foods, which is highlighted in the paper [7], while paper [8], discussed the fundamentals involved in developing a CFD solution giving a state-of-the-art review on various CFD applications in the food industry.The results of numerical simulations are available for different mixing 2D and 3D geometries, the rheological behaviour of fluid and flow regime.Mixing geometries include static mixers [9,10], passive mixer [11], dough kneaders [12][13][14][15], stirred tank reactors [16][17][18][19][20], and continuous mixers [21][22][23].Effects on mixing flows of Newtonian [19,24], and non-Newtonian fluids [24,25] are analysed in all three flow regimes, from laminar [26,27] through transitional [18] to turbulent [28,29].
Power consumption analyses of Newtonian and non-Newtonian fluids in the laminar regime for different mixer configurations were given in papers [3,[30][31][32].
The VOF model of Hirt and Nichols [33], which used finite-difference method was applied in predicting the free surface shape.The effect of rotation speed on mixing performance of a Newtonian fluid in an unbaffled stirred tank with an impeller was investigated by Lamberto et al. [34] using a rotating reference frame.Numerical simulation of mixing in a baffled stirred-tank reactor was performed and presented in paper [35].They started with a multiple reference frame (MRF) based on the steady state algorithm that was followed by full transient sliding mesh approach.The turbulent flow with a free surface vortex in unbaffled vessels agitated by a paddle impeller and a Rushton turbine was analysed in [36].In all these cases it was demonstrated that the free surface was captured well using the standard VOF method, but simulations resulted in smooth, parabolic shaped interfaces with no entertainment of gas into the impeller swept region.The transient mixing processes are also modelled using the mesh superposition technique, a fictitious domain method proposed in [37], or the mesh immersion technique recently proposed by Valette et al. in paper [38] where cumbersome mesh generation at each time step is avoided.
This research was initialized by the manufacturer of the horizontal twin-shaft paddle batch mixer for industrial usage.The aim was to obtain the influence of the angular velocity on mixing of a shear-thinning soy flour-water paste.At the beginning two rotational speeds were defined, with numerical simulations of laminar mixing using the VOF multiphase model and dynamic mesh were conducted.Power consumption, dispersive mixing index, and volume fraction of the air drawn into the paste were calculated.Obtained results were sufficient for the manufacturer to define operational rotational speed, thus additional ones were not considered.

Geometry of the Mixer, Numerical Model and Grid
The geometry of the mixer is shown in Fig. 1a.The initial height of the paste was h p =0.87 m.There were two contra-rotating shafts, each with seven paddles.The origin of the mixer coordinate system is set at the level of the shaft axes in the middle between them.The axes of the shaft were parallel to the x axis with coordinates y= 0 and z = ±340 mm.The mixer was partially filled with the paste so the paddles were not completely immersed in it, while in the rest of the vessel was the air.The multiphase VOF model was used and considerable entertainment of the air into the region with paste was noticed making the mixing process fully unsteady so that the MRF was not able to give reliable results.Also, it was not possible to use sliding meshes due to the mutual position of the paddles and shafts, but instead the dynamic mesh with remeshing of the complete numeric domain at each time step had to be employed.
The geometry of the tank and shafts with paddles was prepared in the AutoCAD software and then exported to GAMBIT as a "step file".After checking and editing the geometry, initial unstructured mesh was generated and exported to the solver.The time averaged shaft torque served for conducting grid independence study, so that the final grid was consisted of 608869 tetrahedral.The mesh in the whole numerical domain is shown in Fig. 1b, and in the vertical plane perpendicular to the x axis, Fig. 1c.

CFD Modelling
Numerical simulations of 3D fully transient flow in the mixer were performed using the commercial computational fluid dynamics (CFD) software ANSYS FLUENT.The volume of fluid (VOF) multiphase model with Eulerian-Eulerian approach was used.The multiphase flow consisted of two phases, the air and non-Newtonian highly viscous paste with 67% of moisture.The VOF formulation is based on the fact that two phases are not interpenetrating and in each cell or control volume, the volume fractions of two phases sum to unity, 1 The interface between the phases is obtained from the solution of a continuity equation for the volume fraction of each of the secondary phase.Herein, Eq. 1 represents the continuity equation written for each of the phases, if the source, S m , was added: Continuing, the momentum equation, shown in Eq. 2, that was solved over the domain, by the VOF model, was tracking the volume fractions of each phase where g ρ is the weight of the rigid body and F is the external force.The energy equation is contained energy shares of phases in the general form of Eq. 3 The energy and analogously the temperature are treated as the mass-averaged values through the VOF application, as represented in Eq. 4 air air air paste paste paste air air paste paste The density and viscosity in each cell are defined by the properties of the phases and their volume fraction as given in Eq. 5.The rheological behaviour of the paste, determined experimentally at 25 0 C, was provided for numerical analysis by the manufacturer.
( 1) (1 ) water water water air water water water air The interface between phases, for the cases of partly filled control volume, 0 1 paste α < < has to be interpolated and for this purpose the geometric reconstruction (Geo Reconstruction) scheme was selected.This scheme leaded the interface reconstruction, adaptation for numerical solving, by linearization of the interface segments, which are predefined by the grid, as shown in Fig. 2. Due to pronounced Weissenberg effect the shear rate has been changed from 0 to 10 s−1.The shear thinning behaviour of the paste obeyed a power law model with a power law index n = 0.2446, a consistency coefficient K = 456 Pa s n , and a coefficient of determination R 2 = 0.997.The paste was modelled as an incompressible fluid with a density equal to 878 kg/m 3 .
No slip conditions were set on the tank walls, shafts and paddles.The open upper surface of the mixer is defined as pressure inlet with zero gauge total pressure (Fig. 1b).The shape of the domain and interface between phases was changing all the time due to the motion of the shafts with paddles.Dynamic mesh was employed in order to obtain the full transient solution.The update of the volume mesh was performed automatically at each time step depending on the new position of the moving boundaries.Two methods were used for the update of the volume mesh, spring-based smoothing and local cell remeshing.The motion of the shafts with paddles was specified as a rigid body motion using userdefined function (UDF) option in the software.
The PISO algorithm was applied for the Pressure-Velocity coupling for transient flow calculations, while the Body Force Weighted (BFW) scheme was used for pressure.The modified HRIC scheme was chosen for volume fraction because it improves the robustness and stability of the calculations and is generally applied in cases with sharp interfaces.The algorithm PISO and schemes BFW and HRIC are included in the ANSYS.Aforementioned algorithm and schemes are recommended in ANSYS FLUENT for this kind of problem.

RESULTS AND DISCUSSION
The unsteady VOF simulations with serial solver were performed for two operating speeds, i.e. 10 and 25 RPM.The implicit first-order accurate time scheme was used by 720 time steps for one revolution.The periodic solution was reached after five revolutions, whereas CPU time was about 280 h with Intel Xeon at 2.66 GHz.

Power consumption
The flow regime is defined by the value of the Reynolds number, defined by Padron (2001) as where ρ is the density of the fluid, N rotational speed and μ viscosity.The average value of the strain rate was about 1, so that high value of the dynamic viscosity was obtained.The largest value of the Reynolds number of the simulations was below 10.The power consumption is usually represented by a power number given as a function of the total power for both shafts, P: The power constant, K P , is defined as The value of the power constant, K P , was about 3.85 based on the average torque.

Dispersive mixing efficiency
The type of flow that occurs in a system is an indication of mixing effectiveness.Generally, three types of flow can appear, from simple shear flow as the most common in food processing, to elongation flow and rotational or transitional flow.An elongation flow that appears when the velocity is steadily increased while the particles are moving along the streamline is difficult to generate and is usually observed only for a short time.Such flows are highly desirable for mixing.In pure rotational or transitional flows there is no deformation of a fluid particle and because of that they are not mixing effectively.The parameter that is usually used in describing the type of flow representing the measure of dispersive mixing effectiveness is a flow number or a mixing index [40] defined as where the rate of magnitudes of strain tensor |D| and vorticity tensor |Ω|, are the symmetric and asymmetric components of the velocity gradient tensor.Flow number gives the relative strength of the pure elongation flow component and ranges from 0 for pure rotation through 0.5 for simple shear to 1.0 for pure elongation.It is known that pure elongation or irrotational flows are more effective for dispersive mixing than flows with a rotational component so that no effective mixing can appear in a system with the purely rotational flow.The mixing index defined in this way is not frame-invariant.However, frame-invariant flow-type indices were found to be computationally difficult to evaluate due to the requirement of higher-mesh densities in the simulation, whereas obtained results were similar to that of the nonframe-invariant index.
The contours of the volume fraction in three vertical planes perpendicular to the shafts are shown in Fig. 3.The volume fraction of the paste and air for 10 RPM is presented at 25%, 50% and 75% of the mixer length (Figs.3a-3c  The goal was a good mixing of the soy paste with the air, which was more favourable for higher angular velocity because more air was drawn into the mixing region.The similar behaviour was also obtained from other cross-sections along the mixer.
In Fig. 6 the portion of the paste with volume fraction equal to 0.9 can be seen.This portion was much larger for 10 RPM than 25 RPM and in both cases it occurred on the shafts and in the gap between paddles and mixer walls.
Wall shear stress contours on the shafts and paddles for chosen angular velocities we can see in Fig. 7.The maximum value appeared at the tips of the paddles and was directly correlated to the angular velocity.
The lower mixing index was obtained at the shafts axes level (y = 0): resulting in poor-mixed zones where flow number decreased even to 0.2.Sudden jumps occurred due to the influence of the neighbouring shaft paddles.Values for y = 200 mm were also high, but somewhat lower than for y = -200 mm, as noticeable in Fig. 8. Herein, Fig. 8 presented the results referred to 10 RPM while in Fig. 8 b to 25 RPM.Generally, the impact of the angular velocity was rather small confirming those results obtained by Cheng and Manas -Zloczower [41].The flow number along three parallel lines in the vertical plane perpendicular to the axes of rotation in the middle of the mixer, i.e. at 50% of the mixer length is shown in Fig. 9 a-c.The lines were at y = -200 mm, y = 0 and y = 200 mm.All values were above 0.5, except in the small region where the rotational flow was dominant and poor-mixing occurred.In this zone close to one shaft that was not swept by the paddles of another shaft at a given time mixing index was reduced to 0.25.Similar curves were obtained in other parallel planes along the mixer.
In this paper level of mixing of the soy flour-water paste with air, power consumption and dispersive mixing index of the soy flour-water paste laminar flow in the part-filled horizontal twin-shaft paddle batch mixer constructed for industrial usage was examined.A CFD fully transient multiphase VOF model with dynamic mesh was used in analysing the mixing process of the non-Newtonian shear thinning fluid for two rotational speeds.

CONCLUSION
For the observed industrial mixer with a horizontal twinshaft paddle batch, 3D transient numerical simulations showed that more air was drawn into the mixing region for higher shaft angular velocity.However, the maximum value of the wall shear stress on the shafts and paddles occurred at the tips of the paddles and was directly correlated to the angular velocity.
The results of the CFD analysis show that the angular velocity changes from 10RPM to 25RPM have not influence on the mixing-index value.
In the active mixing region the high shear flow and mixing index were noticed, but also there were the poor-mixing zones with the almost rotational flow and mixing index as low as 0.2.These zones were continually changing, depending on the mutual positions of both shafts and their paddles.But generally, they were next to the shafts and at both ends of the mixer.The contra-rotating shafts and close clearances between the paddles and between the paddles and barrel walls provide uniform mixing.
A great deal of engineering effort goes into designing and improving mixing processes.It can be concluded that described numerical simulations give useful results for the characterization of local mixing of highly viscous, non-Newtonian fluids in mixers with complex geometry consisted of two contra-rotating shafts with multiple paddles.
The future investigation of industrial mixer with horizontal twin-shaft paddle batch design will be based on the results of this research, including CFD simulations for analysis the influence of the paddle shape and their positions on the shafts, on the volume fraction of the air drawn into the paste, power consumption and flow number.

Figure 1 .
Figure 1.The geometry of the mixer (a): mesh in the whole numerical domain (b): mesh elements in the vertical plane perpendicular to the shafts (c)

Figure 2 .
Figure 2. (a) Real phases' interface shape, (b) Interface reshaped by the linearization over the segments

Figure 3 .Figure 4 .
Figure 3.The volume fraction of the paste and air in the vertical plane) after 5 rotations of the shaft.planes for 10 and for 25 rpm at three positions (25%, 50% and 75% of the mixer length) ).The volume fraction contours for 25 RPM at the same positions are presented in Figs 3d-3f.Figs. 3 show the contours after five full rotations of the shaft and paddles.Total mixing times for 10 RPM and 25 RPM are 30 s and 12 s, respectively.The dynamics of mixing during the first rotation of the shaft at 10 RPM is shown for the various angles in Fig. 4. Contours of volume fraction are presented for the same position as in Fig 3. Summarized isometric images of volume fraction, through horizontal and vertical planes [40], at x=0.25m and Θ= 630° while rotating with 10 RPM after 10,5 s, are presented in Fig.5, respectively.

Figure 5 .Figure 6 .
Figure 5. Summarized isometric images of volume fraction, through horizontal and vertical planes at x=0.25m and Θ= 630° while rotating with 10 RPM after 10,5 s

Figure 7 .
Figure 7. Shear stress on the surface of the shaft and its paddles, at 10 RPM (a) and at 25 RPM (b)