Numerical Modeling of the Mixing of Highly Viscous Polymer Suspensions in Partially Filled Sigma Blade Mixers

This paper presents a non-isothermal, non-Newtonian Computational Fluid Dynamics (CFD) model for the mixing of a highly viscous polymer suspension in a partially filled sigma blade mixer. The model accounts for viscous heating and the free surface of the suspension. The rheological model is found by calibration with experimental temperature measurements. Subsequently, the model is exploited to study the effect of applying heat both before and during mixing on the suspension’s mixing quality. Two mixing indexes are used to evaluate the mixing condition, namely, the Ica Manas-Zlaczower dispersive index and Kramer’s distributive index. Some fluctuations are observed in the predictions of the dispersive mixing index, which could be associated with the free surface of the suspension, thus indicating that this index might not be ideal for partially filled mixers. The Kramer index results are stable and indicate that the particles in the suspension can be well distributed. Interestingly, the results highlight that the speed at which the suspension becomes well distributed is almost independent of applying heat both before and during the process.


Introduction
Attaining and ensuring homogeneity in particle-based slurries is not a trivial task-not least due to the complexity in quantifying homogeneity. The problem only gets more complicated when dealing with highly viscous non-Newtonian fluids. However, applications where homogeneity control is paramount are broad, extending to many industries such as food, pharmaceuticals, concrete and medical devices.
Some of the earliest attempts to obtain a quantification of homogeneity were made by Lacey in 1954 [1] and later by P. V. Danckwerts [2]. Both stated the importance of the sampling procedure, which is essential in determining whether particles in a population were evenly distributed. Other studies have used density as a method for the quantification of homogeneity [3,4]. These approaches have used standard statistics methods such as Coefficient of Variance [3] or Analysis of Variance (ANOVA) [4]. This works if you have appreciable unevenness in the distribution and a significant density discrepancy between components [5].
One of the pioneers within the field of mixing is Ica Manas-Zlaczower [6][7][8][9][10][11]. Her famous dispersive mixing index from 1992 is still used today [12][13][14][15][16][17], along with her distributive Cluster Distribute Index [12,17,18]. There have not been many alternatives to Zlaczower's dispersive mixing index. On the other hand, there are a number of alternatives to the Cluster Distribution Index-for example, the Scale of Segregation, a distributive mixing index that was developed back in 1952 [19] and later applied in Computational Fluid Dynamics (CFD) models by Connelly and Kokini in 2007 [20]. Lastly, there is the Lacey mixing index [21], which has been modified over the years. Other indices that have been developed based on the Lacey mixing index include the Kramer mixing index [22] and the Ashton-Valentin mixing index [23].
Heating is an important aspect of mixing, since it can reduce viscosity drastically [12,17,24,25]. The advantage of reducing the viscosity is that it leads to an increased Reynolds Number, which improves mass transfer (i.e., mixing) [25]. The viscous dissipation of heat is a phenomenon that is often neglected in fluid dynamics [25]. However, when dealing with highly viscous fluids, viscous heating can be present if they are mixed with a certain force. Consequently, this leads to an impact on the energy balance and should therefore only cautiously be ignored, as shown in previous studies [12,26,27].
Sigma blades can often be used as the rotating mixing element in a mixer when dealing with highly viscous fluids [28]. They can be used either as a single-element [3] or, as often seen, as a twin-arm mixer [29]. Sigma blades have previously been simulated by Connelly [30,31] and evaluated with respect to the dispersive mixing index. However, to the best knowledge of the authors, no numerical models have been presented in the literature that account for a free surface, a non-Newtonian material behavior and viscous dissipation when mixing with Sigma blades. This paper presents a new non-isothermal, non-Newtonian CFD model for the sigma blade mixing of an adhesive suspension, where both the free surface and viscous heating are accounted for. The suspension is modeled as a viscoplastic fluid, and model calibration is performed via optical temperature measurements. To evaluate the mixing quality, we have used Zlaczower's dispersive mixing index as well as Kramer's distributive mixing index. The model is used to investigate the effect of applying heat both before and during the process on the mixing quality. The rest of the paper is organized as follows: Section 2 introduces the experimental setup as well as the numerical model, while in Section 3, the results are presented and discussed. Section 4 summarizes the conclusions of the study.

Experiments
The mixing material consists of a highly viscous polymer fluid and a multi-component powder blend. The powder contains both colloidal and microscale particles. The specific fluid/powder suspension is a piece of intellectual property and cannot be disclosed. The viscosity of suspensions changes with the volume fraction of the powder [32,33], but this is not accounted for in the numerical model, as it is assumed to have a limited effect on the results. The physical data, besides viscosity, are shown in Table 1. The density measurements were performed with a Sartorius YDK03 Density Kit for Analytical Balances, which uses the Archimedes principle. The heat capacity and conductivity measurements were performed on a TCi-3-A from C-Therm Technologies Ltd., which uses the Modified Transient Plane Source (MTPS) method that has been used in textile research [34]. Heat capacity 1389 The experimental mixing was carried out on a Linden LK II 1 machine; see Figure 1. The machine has a mixing chamber with two sigma blades rotating in opposite directions. The front blade rotates at 60 rpm, and the one behind rotates at 19 rpm. Historically, this has been a typical setting for highly viscous mixing [28]. The mixer also has a control panel to switch the rotation direction. In addition, there is a digital display that indicates the Polymers 2023, 15,1938 3 of 12 temperature, as measured at the bottom of the mixer. Attached to the mixing equipment is a heat exchanger that adjusts the temperature for the mixing process. A vacuum pump is also mounted to ensure that air is not trapped inside the mixture. Two walls in the mixer can provide heat to the system. If this function is on, the temperature will be 80 • C. The experimental mixing was carried out on a Linden LK II 1 machine; see Figure 1. The machine has a mixing chamber with two sigma blades rotating in opposite directions. The front blade rotates at 60 rpm, and the one behind rotates at 19 rpm. Historically, this has been a typical setting for highly viscous mixing [28]. The mixer also has a control panel to switch the rotation direction. In addition, there is a digital display that indicates the temperature, as measured at the bottom of the mixer. Attached to the mixing equipment is a heat exchanger that adjusts the temperature for the mixing process. A vacuum pump is also mounted to ensure that air is not trapped inside the mixture. Two walls in the mixer can provide heat to the system. If this function is on, the temperature will be 80 °C.

Numerical Model
The CFD model simulates the mixing of the suspension in the mixing chamber and is developed in the commercial software FLOW-3D, which has successfully simulated other processes such as casting [35,36] and 3D printing [37,38] highly viscous fluids. An illustration of the mixer geometry is seen in Figure 2. At t 0, the fluid is well distributed at the bottom, and the powder is placed on top of the fluid. At t 0, the sigma blades start rotating. A no-slip boundary condition is applied on all solid surfaces. The computational domain is meshed with a uniform grid consisting of ~50.000 elements, which was arrived at after a mesh sensitivity analysis. The two side walls perpendicular to the sigma blades apply heat to the system; see Figure 2.

Numerical Model
The CFD model simulates the mixing of the suspension in the mixing chamber and is developed in the commercial software FLOW-3D, which has successfully simulated other processes such as casting [35,36] and 3D printing [37,38] highly viscous fluids. An illustration of the mixer geometry is seen in Figure 2. At t = 0, the fluid is well distributed at the bottom, and the powder is placed on top of the fluid. At t > 0, the sigma blades start rotating. A no-slip boundary condition is applied on all solid surfaces. The computational domain is meshed with a uniform grid consisting of ∼50,000 elements, which was arrived at after a mesh sensitivity analysis. The two side walls perpendicular to the sigma blades apply heat to the system; see Figure 2. The flow is transient and non-isothermal since the viscosity is temperature dependent. The material is modeled as an incompressible substance, and thus, the density is approximated as constant. Hence, the flow is computed by considering the mass, momentum and energy conservation: The flow is transient and non-isothermal since the viscosity is temperature dependent. The material is modeled as an incompressible substance, and thus, the density is approxi- where ρ is the density, p is the pressure, k is the thermal conductivity C p , is the specific heat capacity and q is the heat flux vector. τ is the material deviatoric stress tensor and is calculated as τ = 2µ . γ, T D. D is the trace of the deformation rate tensor, and it is . γ is the shear rate and is calculated by The gravitational acceleration, G, is given by (0, 0, −9.82) m s 2 . The software used the finite volume method to discretize the governing equations. The equation of energy is calculated explicitly, while the viscous stress and pressure are solved implicitly. The advection is solved explicitly with first-order accuracy. The free surface is calculated with the volume of fluid technique [39].
The material behaves as a non-isothermal viscoplastic fluid and is simulated by a modified Carreu model: where µ 0 is the zero-share-rate viscosity, µ ∞ is the infinity-share-rate viscosity, n is the power exponent, λ is the time constant, E is an energy function which is dependent on the fluid temperature, a and b are empirical constants and T re f is the reference temperature.
The applied values of the variables in Equations (5) and (6) are seen in Table 2. The shear rate-dependent variables are obtained via a rheological characterization of the fluid; see Figure 3. The temperature-dependent variables are obtained by a calibration, as the fluid was too viscous to perform measurements at low temperatures. The calibration is reported in Section 3.1. Table 2. Viscosity data of the simulated fluid.

Symbol
Value Unit  The evaluation of the mixing is carried out by a dispersive-and distributive-mixing index. The dispersive mixing is quantified through the Manas-Zlacower mixing index, λ ; see Equation (8).
is the vorticity. When λ is close to 0, the mixing is purely rotational driven, while at 0.5 and 1, the mixing is shear-driven and purely elongationdriven, respectively. The latter leads to a desirable faster breakup of agglomerates [11,13,29].
The second approach is a distributive mixing index where an artificial concentration or material is used. The concentration is defined as a scalar with zero diffusion, which does not affect the physics (such as the viscosity and density). Thus, the only effect is pure mixing. The mean of the concentration, ̅ , is 0.5. The dimensionless concentration, , and the dimensionless variance across the domain that contains fluid, , are found by Equations (8) and (9). From the variance, the Kramer dispersive mixing index, , can be found by Equation (10)   The evaluation of the mixing is carried out by a dispersive-and distributive-mixing index. The dispersive mixing is quantified through the Manas-Zlacower mixing index, λ mz ; see Equation (8). ω is the vorticity. When λ mz is close to 0, the mixing is purely rotational driven, while at 0.5 and 1, the mixing is shear-driven and purely elongationdriven, respectively. The latter leads to a desirable faster breakup of agglomerates [11,13,29].
The second approach is a distributive mixing index where an artificial concentration or material is used. The concentration is defined as a scalar with zero diffusion, which does not affect the physics (such as the viscosity and density). Thus, the only effect is pure mixing. The mean of the concentration, c, is 0.5. The dimensionless concentration, c i , and the dimensionless variance across the domain that contains fluid, S 2 , are found by Equations (8) and (9). From the variance, the Kramer dispersive mixing index, M Kramer , can be found by Equation (10) Table 3 presents the process parameters for the simulations that are studied in this paper.

Calibration
In the experimental setup, an optical thermal probe was placed inside the mixing chamber to measure the temperature of the suspension. The probe logged a surface temperature every 10 s. The initial temperature of the fluid was measured to be 23 • C, and no heat was applied while mixing. The experimental findings and the results of the calibration simulation, CA23, are shown in Figure 4. The experimental measurements show that viscous heating is present and that the temperature increases most rapidly during the first 250 s. The fluctuation around 250 and 500 s is possibly due to the brushed steel occasionally being measured instead of the fluid. The parameters a and b in Table 2 are calibrated in order to make the CFD model predict the same temperature evolution as that seen in the experiment. By doing so, the experimental findings and the numerical results are in very good agreement. The absolute mean error between the two is 1.72 • C. The rheological model obtained by calibration is presented in Figure 5 for the temperatures 23, 50 and 80 • C.

Temperature Distribution
The temperature distribution is close to uniform when mixed under adiabatic boundary conditions, but this is not the case when applying the non-adiabatic boundary condition, as seen in Figure 6, which illustrates the temperature profile for NA50 at 15 s. The profile shows that the temperature is highest near the heated walls. The temperature decrease along the x-axis is due to the fluid not being exposed to the fixed wall temperature. The average temperature increases with time due to the viscous heating and heated walls.

Temperature Distribution
The temperature distribution is close to uniform when mixed under adiabatic boundary conditions, but this is not the case when applying the non-adiabatic boundary condition, as seen in Figure 6, which illustrates the temperature profile for NA50 at t = 15 s. The profile shows that the temperature is highest near the heated walls. The temperature decrease along the x-axis is due to the fluid not being exposed to the fixed wall temperature. The average temperature increases with time due to the viscous heating and heated walls.  23. The simulations with applied wall temperature increase faster in terms of temperature as compared to their counterpart with adiabatic boundary conditions, except for NA80, where the walls effectively cool the suspension. A23 and A50 do not have enough time to reach 80 °C, while A80 reaches a temperature slightly higher than 80°.  Figure 7 for all simulations except CA23. The simulations with applied wall temperature increase faster in terms of temperature as compared to their counterpart with adiabatic boundary conditions, except for NA80, where the walls effectively cool the suspension. A23 and A50 do not have enough time to reach 80 • C, while A80 reaches a temperature slightly higher than 80 • C. The average temperature of the fluid as a function of time is illustrated in Figure 7 for all simulations except 23. The simulations with applied wall temperature increase faster in terms of temperature as compared to their counterpart with adiabatic boundary conditions, except for NA80, where the walls effectively cool the suspension. A23 and A50 do not have enough time to reach 80 °C, while A80 reaches a temperature slightly higher than 80°.

Velocity Field
The velocity field gives a good indication of how the flow pattern is inside of the mixer. In Figure 8, the velocity field for A80 is presented at different times. As expected, the sigma blade at the lowest x-value is the one with the highest velocity, while the second sigma blade spins slower to ensure mixing in the whole system. Noticeably, some of the fluid sticks to the wall, which is due to the viscoplastic effect of the fluid and the no-slip condition on the walls.

Velocity Field
The velocity field gives a good indication of how the flow pattern is inside of the mixer. In Figure 8, the velocity field for A80 is presented at different times. As expected, the sigma blade at the lowest -value is the one with the highest velocity, while the second sigma blade spins slower to ensure mixing in the whole system. Noticeably, some of the fluid sticks to the wall, which is due to the viscoplastic effect of the fluid and the no-slip condition on the walls.

Dispersive Mixing
In Figure 9, the histogram results of at 1800 s are shown for all simulations except 23. The results are fairly similar for all simulations, which indicates that preheating and heating during the process do not have a big influence on the mixing quality after 30 min. The dispersive mixing index is primarily between 0 and 0.5, which is similar to the findings of Ahmed and Chandy [12], who simulated a two-wing rotors mixer. However, the values are lower than the predictions by Connelly and Kokini [30], who

Dispersive Mixing
In Figure 9, the histogram results of λ mz at t = 1800 s are shown for all simulations except CA23. The results are fairly similar for all simulations, which indicates that preheating and heating during the process do not have a big influence on the mixing quality after 30 min. The dispersive mixing index is primarily between 0 and 0.5, which is similar to the findings of Ahmed and Chandy [12], who simulated a two-wing rotors mixer. However, the λ mz values are lower than the predictions by Connelly and Kokini [30], who studied a fully filled sigma blade mixer, which indicates that a partially filled mixer seems to have a decreasing effect on how well dispersed the particles will be in the suspension. Note that the presence of λ mz values between 0 and 0.1 is mainly due to the viscoplastic effect that makes the suspension more prone to have zones with a shear rate of zero and therefore limited mixing.  Figure 10 shows the mean value, , illustrated at different times. It is seen that the mean value fluctuates by ~10%. The fluctuations are most likely due to the free surface of the mixture (i.e., the mixer is partially filled).

Distributive Mixing
In Figure 11, the Kramer mixing index is illustrated at different time values for all simulations except CA23. The results show that all simulations obtain a high value (i.e., close to 1) after 2 min, which indicates that the particles in the suspension end up being well distributed. A80 and NA80 reach a high value the fastest (due to the lower viscosity at high temperatures), but this is marginally faster than the rest, thereby highlighting that energy potentially can be saved by not preheating the mixture Figure 9. Histogram of the λ mz at t = 1800 s. Figure 10 shows the mean value, λ mz , illustrated at different times. It is seen that the mean value fluctuates by~10%. The fluctuations are most likely due to the free surface of the mixture (i.e., the mixer is partially filled).  Figure 10 shows the mean value, , illustrated at different times. It is seen that the mean value fluctuates by ~10%. The fluctuations are most likely due to the free surface o the mixture (i.e., the mixer is partially filled).

Distributive Mixing
In Figure 11, the Kramer mixing index is illustrated at different time values for al simulations except CA23. The results show that all simulations obtain a high value (i.e., close to 1) after 2 min, which indicates that the particles in the suspension end

Distributive Mixing
In Figure 11, the Kramer mixing index is illustrated at different time values for all simulations except CA23. The results show that all simulations obtain a high M Kramer value (i.e., close to 1) after 2 min, which indicates that the particles in the suspension end up being well distributed. A80 and NA80 reach a high M Kramer value the fastest (due to the lower viscosity at high temperatures), but this is marginally faster than the rest, thereby highlighting that energy potentially can be saved by not preheating the mixture without compromising the mixing of the suspension. Similarly, the results show that no substantial gain will be obtained regarding the mixing time by applying heat during the process. One should keep in mind that at a fixed mixing velocity, the Sigma blades require more energy at low temperatures due to the higher viscosity of the material. In addition, the Sigma blades will be exposed to larger forces at high viscosities, which can affect the blades' lifespan. without compromising the mixing of the suspension. Similarly, the results show that n substantial gain will be obtained regarding the mixing time by applying heat during th process. One should keep in mind that at a fixed mixing velocity, the Sigma blades requir more energy at low temperatures due to the higher viscosity of the material. In addition the Sigma blades will be exposed to larger forces at high viscosities, which can affect th blades' lifespan.

Conclusions
This work presents a non-isothermal, non-Newtonian CFD model that simulates partially filled Sigma blade mixer by taking into account viscous heating and the free sur face of the suspension. The model is calibrated against physical temperature measure ments, and good agreement between the experimental and numerical results is obtained The results also highlight the importance of accounting for viscous heating when mode ing the mixing of the suspension at hand. A numerical parameter study (i.e., varying wa temperature and pre-heating) was made in order to improve the understanding of th mixing conditions. The results of the model showed that the dispersive mixing index, was primarily between 0 and 0.5, and it was not affected much by changing the wall an pre-heating temperature. Fluctuations of approximately 10% were found in the dispersiv mixing index for most simulations, which was associated with the partially filled mixe (i.e., free surface of the suspension). This can be seen as a limitation of the mixing index The model predicted via the Kramer mixing index that the suspension in all simulation would be well distributed after 1800 s. In addition, the results illustrated that pre-heatin as well as heating during the process did not provide a substantial gain when it came t how fast a high Kramer index would be obtained. This illustrated that one could poten tially save energy and time by eliminating the heating. This was a very interesting finding as, intuitively, one could have expected that heating would lead to lower viscosity an

Conclusions
This work presents a non-isothermal, non-Newtonian CFD model that simulates a partially filled Sigma blade mixer by taking into account viscous heating and the free surface of the suspension. The model is calibrated against physical temperature measurements, and good agreement between the experimental and numerical results is obtained. The results also highlight the importance of accounting for viscous heating when modeling the mixing of the suspension at hand. A numerical parameter study (i.e., varying wall temperature and pre-heating) was made in order to improve the understanding of the mixing conditions. The results of the model showed that the dispersive mixing index, λ mz , was primarily between 0 and 0.5, and it was not affected much by changing the wall and pre-heating temperature. Fluctuations of approximately 10% were found in the dispersive mixing index for most simulations, which was associated with the partially filled mixer (i.e., free surface of the suspension). This can be seen as a limitation of the mixing index. The model predicted via the Kramer mixing index that the suspension in all simulations would be well distributed after 1800 s. In addition, the results illustrated that pre-heating as well as heating during the process did not provide a substantial gain when it came to how fast a high Kramer index would be obtained. This illustrated that one could potentially save energy and time by eliminating the heating. This was a very interesting finding, as, intuitively, one could have expected that heating would lead to lower viscosity and therefore substantially faster mixing. In future work, focus will be put on extending the model to account for the effect of powder concentration variations on the viscosity and thermal conductivity, as this can especially affect the initial phase of the mixing.