Abstract
Rubber mixing is an important link in the production of rubber products. Computational fluid dynamics (CFD) simulation is often used to explore the effect of rubber mixing parameters on rubber mixing effect. Previous CFD-based rubber mixing simulation studies did not consider the impact of using 2D or 3D numerical calculation models on the numerical simulation results. In order to investigate the differences between 2D and 3D numerical computational models in rubber compounding CFD simulation problems, in this paper, we compare and analyze the results obtained from 2D and 3D computational models under different rotational speed conditions to investigate the differences between the models in the numerical simulation of rubber compounding. Three different experimental speeds of the rubber mixer—39, 44, and 49 r/min—were set during the study using 2D and 3D asynchronous rotor models with a speed ratio of 1.15, respectively. The rubber was processed using the Bird–Carreau model. The phase interface between rubber and air was calculated using the volume of fluid (VOF) method. The numerical simulation results of different models show that the rotational speed set to 49 r/min shows the best dispersion distribution effect; the mixing effect and speed change rule obtained by the 2D model are consistent with the results obtained by the 3D model. The performance of the results of the two models is consistent when exploring the numerical simulation of rubber compounding.
1. Introduction
Rubber products are widely used in the automotive industry, such as automobile tires, window seals, automobile suspensions, and so on. Before rubber is processed into various products, the raw material has to go through many procedures, and rubber mixing is an important part of the rubber product processing. Rubber mixing process parameters, such as rubber filling coefficient and rubber mixer speed, will directly affect the final product quality of rubber products [1]. Rubber mixing process parameters are usually determined using the test method, but the test method has the disadvantages of a long cycle time, high cost, and low efficiency. In recent years, with the development of supercomputers and the proposal of more accurate numerical simulation algorithms, the efficient advantages of numerical simulation analysis have become more and more significant. When using numerical simulation methods to study rubber mixing problems, the computational models can be categorized into 2D planar models and 3D spatial models. The 2D model has only two dimensions, and the Navier–Stokes equations have fewer unknown variables, resulting in a faster solution. Although the 3D model matches the actual situation, the 3D model has three dimensions, and the number of unknown variables in the Navier–Stokes equation increases, resulting in a slower solution speed. Therefore, it is necessary to explore the differences between 2D and 3D computational models in numerical simulation calculations of rubber compounding and to clarify the influence of the two models on the calculation results.
At present, some scholars have conducted some studies on the numerical simulation process of rubber mixing. Song et al. [2] used a 2D model to study the effect of different rotor structures on the mixing effect of rubber. Das et al. [3] investigated the effect of different speed ratios on the rubber mixing effect under isothermal conditions using 3D modeling and indicated that the optimum speed ratio was 1.125. Connelly and Kokini [4] used 3D modeling to study the rotor mixing effects of Sigma structures. Kim and White [5] experimentally investigated the optimal filling factor of different rubbers, and the experiment showed that the rubber mixing effect was better when the filling factor was 75%. Frungieri et al. [6] used a combined CFD and DEM approach to study the dispersion process of rubber agglomerates in a 2D model. Han and He [7] used a 3D model to investigate the difference between partial and complete filling under isothermal conditions. Liu et al. [8] analyzed the transient distribution position of rubber in a mixer using a 2D model and obtained results similar to the conclusions drawn by Han and He. Poudyal et al. [9] analyzed the rubber mixing effect of isothermal and nonisothermal conditions using a 2D model and concluded that the change rule of the evaluation indexes under the two conditions is consistent. Although related scholars have explored the numerical simulation process of rubber mixing, the study was carried out on the basis of selected 2D or 3D computational models and did not consider whether the difference in model selection would have a greater impact on the simulation calculation results.
When performing numerical simulation calculations, although the analysis using 3D models will be more in line with the actual situation, the whole calculation process will consume a lot of time, which is not conducive to the exploration of the relevant laws. In contrast, the calculation efficiency of a 2D model is higher than that of a 3D model. Therefore, some scholars tend to use 2D models for analyzing and calculating [10, 11]. However, some scholars tend to use the 3D model for analyzing and calculating because it is closer to the actual working conditions [12, 13]. Different scholars used different models to explore the research. First, it was determined to use 2D or 3D computational models, according to which the effects of rotational speed, speed ratio, filling factor, and phase angle on the rubber mixing effect were investigated, as shown in Figure 1. Obviously, this process does not consider the effect of model differences on the results, and the study of the effect of 2D and 3D models on the final rubber mixing effect is rarely reported. Therefore, the purpose of this paper is to investigate the effect of different models on the mixing effect by simulating the rubber mixing process of the two computational models at different rotational speeds and evaluating the rubber mixing effect using the Lagrangian statistical method [14], to investigate the consistency of the conclusions obtained from 2D and 3D computational models at different rotational speeds, and to provide theoretical and methodological guidance for improving the quality of rubber compounding.
2. Geometry and Materials
The geometric model in this paper consists of a mixing chamber and two biplane nonmeshing rotors. Figure 2(a) shows the complete 3D rotor model, and Figure 2(b) shows the cross-section of the 3D rotor model at Z = 473 mm. Figure 3 shows the 2D rotor model. The starting position of the rotor is the same for both the 2D rotor model and the 3D rotor model for the Z = 473 mm section. The left and right rotors rotate in opposite directions, with the left rotor rotating clockwise and the right rotor rotating counterclockwise. The left and right rotors are set as asynchronous rotors.
(a)
(b)
In the mixing chamber, set the rest of the area to be air, except for the rubber material. Air has a density of 1.225 kg/m3 and a viscosity of 1.78 × 10−5 kg/(m s). The material parameters of rubber were measured experimentally, and the rubber material parameters are shown in Table 1. In order to characterize the flow properties of rubber in a rubber mixer, the rheological properties of rubber materials are characterized using the Bird–Carreau model, which behaves as a pseudoplastic fluid at high shear rates and as a Newtonian fluid at low shear rates [15]. The relationship between shear rate and viscosity can be expressed as follows:where is zero shear stress, is infinite shear stress, is time constant, and n is power–law index.
3. Basic Equation
For polymer melt flow, the basic governing equation can be expressed as follows:
Continuity conservation equation:where is density of the fluid; is time; and are components of the velocity vector in the x, y, and z directions.
Momentum conservation equation:
Momentum equation in the X direction:
Momentum equation in the Y direction:
Momentum equation in the Z direction:where is pressure of the fluid; is the cell surface stress ; and are the mass force along the x, y, and z directions on the surface of the unit.
Since the research in this paper is carried out under isothermal conditions, the energy equation is not considered.
On the other hand, the simulation in this paper is partially filled, and the mixing chamber is composed of rubber and air, so it is necessary to track the rubber–air phase interface. Here, phase interface tracking is performed using the Euler-based multiphase flow volume of fluid (VOF) model. The two different phases (air and rubber) are controlled by a set of continuity and momentum equations, but the volume fraction of each phase can be tracked in every cell over the entire region. The equations of the VOF model are as follows:where is volume fraction of phase m.
In order to analyze the rubber mixing effect, a group of massless particles is set, and the tensile length, distribution index, and cumulative probability distribution of the maximum shear stress are calculated based on the massless particles. Track the position of particles during a numerical simulation. In order to track these particles, a velocity field interpolation scheme is used here, defined as follows:where is velocity of the particle.
4. Calculation Strategy
The 2D and 3D computational models were drawn using CATIA software, and the computational models were preprocessed through ANSYS DesignModeler software. The finite element model meshing is done using Fluent meshing. Through the 2D and 3D model mesh independence tests and considering the computational efficiency, the mesh is divided in the following way: For the static region of the 3D model, it is divided into a mixed mesh of 7 mm tetrahedra and hexahedra, with a total of 159,473 mesh cells. For the dynamic region, a polyhedral mesh of 6–7 mm is used with 1013,137 meshes, in which the boundary layer mesh of the rotor boundary region is four layers and the growth rate of the boundary layer is 1.2. The static region is divided into a hybrid mesh of hexahedra and tetrahedra with 146,336 meshes, as shown in Figure 4. The static region of the 2D model uses a 4-mm triangular grid with a total of 6,934 meshes; the dynamic region uses a hybrid triangular quadrilateral grid with encrypted meshes at the tip of the rotor; the 2D model has a total of 21,700 grid cells, as shown in Figure 5.
Numerical simulation calculations were performed using the commercial CFD software Fluent with partially filled isothermal conditions. The sliding grid method is used to realize the motion process of the rotor. The dynamic and static regions exchange data through the grid intersection interface. Since rubber is a highly viscous fluid and its flow process is laminar, the flow model is chosen as a laminar flow model [16]. The pressure–velocity coupling is implemented using the PISO algorithm, momentum is discretized using the second-order format, pressure is discretized using the PRESTO! algorithm, and the time term is discretized using the first-order implicit format, with the time step set to 0.0005 s.
The volume fraction of the rubber phase was set to 75%, and the volume fraction of the air phase was set to 25%. The numerical simulation experimental speeds were set to 39, 44, and 49 r/min for the left rotor and 44.85, 50.6 , and 56.35 r/min for the right rotor. Figure 6 shows the initial distribution of the rubber inside the rubber mixer.
5. Results and Discussion
After the numerical simulations, the dispersion mixing effect and distribution mixing effect of the rubber are analyzed using the Lagrangian statistical method. Here, dispersion mixing denotes the process of breaking large agglomerates of rubber polymer into small particles. Distributed mixing is the process of uniform distribution of additives and dispersed small particles of carbon black throughout the rubber polymer [17]. Existing scholars have pointed out that the dispersive mixing is realized by the shear and tensile forces generated by the rotary motion of the rotor, which is quantitatively analyzed using two indexes: the mixing index and the cumulative probability distribution of the maximum shear stress [18]. The effect of distributed mixing was quantitatively analyzed using two metrics: the distribution index and the mean length of stretch [19]. By comparing the pressure distributions, rubber phase distributions, mixing indices, cumulative probability distributions of maximum shear stresses, distribution indices, and average tensile lengths calculated by 2D and 3D computational models at different rotational speeds, we explored the differences between 2D and 3D computational models in the numerical simulation of the rubber mixing problem, as well as the effects of the differences in the model choices on the results of numerical calculations.
5.1. Distribution of Pressure Field
Figure 7 shows the distribution of pressure areas for the 3D model and 2D model at different rotational speed conditions. The region in front of the rotor rotation direction belongs to the high-pressure gradient region, the region behind the rotor belongs to the negative pressure gradient region, and the maximum pressure region is concentrated in the gap region formed by the rotor tip and the wall of the mixing chamber. This is consistent with the results of Han [7]. It shows the effectiveness of the numerical simulation method in this paper. During the rubber mixing process, the rubber flows into the channel formed between the outer surface of the rotor and the wall of the mixing chamber. The rubber flows from the front area of the rotor into the high shear stress zone formed by the top of the rotor and the wall of the mixing chamber and then out to the rear area of the rotor. During this process, the large rubber aggregates are broken up and dispersed into smaller aggregates. After analyzing the percentage of high-pressure regions in the 2D and 3D models, as shown in Table 2, it was concluded that a large pressure gradient occurs when the rotational speed of both models is set to 49 r/min. This stratified pressure gradient facilitates the flow of rubber in the high shear stress region and improves the crushing and dispersing effects of large-grained rubber aggregates. Scholar Das et al. [20] found in his study that increasing the rotational speed leads to an increase in the pressure value, which is in agreement with the results obtained in this paper and proves the validity of the numerical simulation method in this paper. Comparing the pressure area distribution between the 2D and 3D models, the pressure in the mixing chamber of the 3D model is higher than the pressure field in the mixing chamber of the 2D model. The reason for this phenomenon is the presence of axial flow in the 3D model, which exacerbates the squeezing action between the rubber materials and produces a higher-pressure field.
Figure 8 shows the trend of the centerline pressure values of the 2D model and the centerline pressure values of the 3D model. The trend of the 2D model and the 3D model is the same. However, at a rotational speed of 44 r/min, the distribution of centerline pressure values for the 2D model and 3D model is opposite to the distribution pattern of pressure values for the remaining two rotational speed conditions. The reason for this phenomenon is that the 3D model only calculates the pressure distribution on the cross-section at Z = 473 mm, and the global pressure distribution of the 3D model cannot be reflected by this cross-section alone.
(a)
(b)
(c)
5.2. Rubber Phase Volume Fraction
Figure 9 shows the cloud view of the rubber phase distribution at 24 s of simulation. A volume fraction of 1 represents a pure rubber phase, and a volume fraction of 0 represents a pure air phase. The volume fraction can also be used as an index for a preliminary evaluation of the mixing effect. When the rotational speed was set to 39 r/min, a large number of air phase agglomerates existed in the mixing chamber. When the rotational speed was set to 49 r/min, the air phase agglomerates decreased, indicating a better rubber mixing effect. In order to be more accurate, by comparing the average volume fraction of the rubber phase of the 2D model and the 3D model, as shown in Table 3, it was concluded that the average volume fraction of the rubber phase was the highest when the rotor speed was 49 r/min. This indicates that the rubber is more fully mixed internally at a speed of 49 r/min. However, it should be emphasized that this analysis is only from a visual point of view and is not a basis for judging the final mixing results.
5.3. Dispersion Mixing
During the rubber mixing process, large-grained rubber polymers are broken into small-grained polymers by shear and tensile stresses. It has been shown that tensile stress can promote the mixing of rubber agglomerates more effectively than simple shear stress [21]. The mixing index is an important indicator to quantify the dispersion effect, and the mixing index is calculated as follows [22]:where |γ| is tensor of deformation rate and |ω| is speed rotation.
Mixing index values range from 0 to 1.0 for pure rotational motion, 0.5 for pure shear motion, and 1 for pure tensile flow.
Figures 10 and 11 show histograms of the mixing index probabilities for the 2D and 3D models, respectively, when the simulation is computed up to 24 s. The highest percentage of tensile flow was observed when the 2D model speed was set to 49 r/min. Compared to the remaining two rotational speeds, the proportion of shear and tensile flow increased. Analyzing the mixing index distribution of the 2D model at three simulation speeds, the proportion of mixing index between 0 and 0.45 is relatively small when the speed is set to 49 r/min. The proportion of mixing index between 0.55 and 0.95 is higher. In summary, it can be concluded that when the rotational speed is set to 49 r/min, the dispersion mixing effect is better than the other two rotational speeds. Comparing the mixing index probability histograms at 24 s for the 3D model, the 3D model has the smallest percentage of rotational flow at a speed setting of 49 r/min and a higher percentage of tensile flow components than at a speed setting of 44 r/min. The mixing index distributions for the two different speed conditions are very close to each other when compared to the speed setting of 39 r/min, but the percentage of tensile flow at the speed of 49 r/min is higher than the percentage of tensile flow at the speed of 39 r/min.
Comparing the variation rules of mixing index between 2D and 3D models under different rotational speed conditions, the trend of mixing index is different between 2D and 3D models when the mixing index is in the range of 0.05–0.45, which is mainly due to the lack of axial flow and the enhanced tendency of rotational flow of the rubber in the 2D model. In the region of 0.55–0.95, which is more favorable for rubber dispersion, the trends of mixing index changes are consistent between the 2D and 3D models. It shows better dispersion mixing at 49 r/min.
The crushing of large-grained polymers in rubber relies mainly on the action of high shear stress at the top of the rotor. The maximum shear stress cumulative probability distribution can visualize the rubber macroparticle polymer experiencing the high shear stress region [23]. To evaluate the maximum shear stress cumulative probability distribution, 2,500 massless particles (2D model) and 8,000 massless particles (3D model) were injected at the beginning of mixing. Massless particles are injected into the mixing chamber in the form of circular surfaces, as shown in Figure 12. Figure 13 shows the distribution of massless particles after 1, 4, and 12 s of numerical simulation calculations for different velocity conditions of the 3D model. After the numerical simulation is carried out for 4 s, the distribution of particles is not uniform for the three velocity cases. As the rotor moves, the particles are carried by the rotor tip to various parts of the mixing chamber, and the distribution gradually becomes uniform. The dispersion mixing effect was quantitatively analyzed using the maximum shear stress cumulative probability.
Figures 14 and 15 show the cumulative probability distribution of the maximum shear stress for the 2D and 3D models, respectively. The 2D model assumes that the shear stress of 200 kPa is the reference value, and the speed is set to 49 r/min. Then, about 70% of the particles are subjected to the maximum shear stress of less than 200 kPa, and 30% of the particles are subjected to the maximum shear stress of more than 200 kPa. When the rotational speed is set to 39 r/min, the curve is closest to the upper-left corner, which indicates that the proportion of particles experiencing higher shear stresses is relatively small, which is not conducive to the breaking of large rubber granules. This is not favorable to the crushing of large rubber agglomerates. The maximum shear stress cumulative probability distribution curve should be shifted to the lower right corner as much as possible, which indicates that more particles will experience higher shear stresses, which is beneficial to the dispersion effect. The curve is closer to the lower right corner when the rotational speed of the 3D model is set to 49 r/min, which indicates that the proportion of particles experiencing higher shear stresses is increased, which is more conducive to the improvement of the rubber mixing effect. This is because with the increase in rotational speed, the shear rate generated by the rotation of the rubber mixer rotor will also increase, and the shear stress will become larger, thus improving the mixing effect of rubber. According to the average shear stress change curve of rubber, as shown in Figure 16, Zhang et al. [24] pointed out that the shear stress suffered by rubber is positively correlated with the rotational speed, which proves the effectiveness of the numerical simulation method in this paper.
(a)
(b)
The combined analysis of the mixing index and the cumulative probability distribution of the maximum shear stress of the two models shows that the best dispersion effect of rubber mixing is achieved when the rotational speed is set at 49 r/min. This indicates that increasing the rotational speed is beneficial to improving the dispersion effect of rubber mixing. This indicates that increasing the rotational speed is beneficial to improving the dispersion effect of rubber mixing. The 2D and 3D models are consistent in investigating the effect of rubber mixing speed on the mixing effect. Scholar Thongpin et al. [25] experimentally investigated the effect of rotational speed on the mixing effect and pointed out that the size of the dispersed phase decreased with the increase in rotational speed, indicating that the dispersion effect became better. Scholar Andreas [26] pointed out that the dispersion effect increases with the increase in rotational speed for the same filling factor, and the dispersion effect decreases after reaching the optimal rotational speed. Once again, the effectiveness of the numerical simulation method in this paper is proven.
5.4. Distribution Mixing
Determining whether the distribution of small particle agglomerates is uniform or not is one of the methods to evaluate the mixing effect. In order to judge the mixing effect of the rubber, the massless particles in the 2D and 3D models were tracked during the numerical calculations. Figure 17 shows the position distribution of the massless particles after 24 s of simulation calculations at different experimental speeds.
Figure 17 shows that the massless particles have been dispersed essentially homogeneously throughout the mixing chamber. In order to make a quantitative judgment of the rubber mixing effect, the distribution index and mean length of stretch were used for quantitative analysis. The distribution index characterizes the extent to which the actual particle position distribution calculated by numerical simulation deviates from the ideal distribution and is calculated by:where is total number of particles, is correlation coefficient between the particle pairs ranging from distance to , and is probability density function such that area under the curve of is always constant and equal to 1.
The distribution index represents the normalized residual value of the deviation from the ideal distribution, and the smaller the value of the distribution index, the closer the actual distribution is to the ideal distribution curve.
Figure 17 shows the actual particle-pair distance distribution curves and the ideal particle-pair distance distribution curves for the 2D model at different moments under different experimental rotational speed conditions. The particle distance distribution curve deviates from the ideal curve at the beginning of mixing, and with the increase in time, the actual particle pair distance distribution curve is infinitely close to the ideal distribution curve.
Figures 18 and 19 show the curves of the distribution indices of the 2D and 3D models with time, respectively. The distribution index is the highest at the beginning of the simulation, and with the change of time, the distribution index shows a decreasing trend, indicating that the actual particle distribution state is getting closer and closer to the particle distribution state in the ideal state. The curves show oscillations caused by the circulation of particle clusters in the mixing chamber. The overall distribution index of the 2D model is at a lower value when the rotational speed is set to 49 r/min, indicating that the actual particle distance distribution curve is closer to the ideal particle distribution curve. It has a better dispersing and mixing effect. Comparing the change process of distribution index of 3D model, the trend curve of distribution index of 3D model is similar to that of 2D model, which indicates that 2D and 3D models have consistency in exploring the influence of rubber mixing speed on mixing effect.
In addition to the distribution index, the distribution mixing effect can also be evaluated using the mean stretch length. The mean stretch length is defined as follows and is used to represent the ratio of different particles to the mean distance [27].where is average distance between particle pairs at moment t and is average distance between particle pairs at the initial moment.
As the mixing process proceeds, the particles begin to move, the distance between pairs of particles changes, and the tensile length increases. The greater the value of the average stretch length, the more violent the movement of the particles, which favors the mixing effect of the rubber.
Figure 20 shows the course of the mean length of stretch of the particles with time. Due to the oscillation of the curve, it is difficult to reflect the effect of rotational speed on the mixing effect. Therefore, we used the cumulative average over a period of time to calculate the average stretching length after filtration. Figures 21 and 22 show the average tensile length after filtration at different rotational speeds for the 2D and 3D models, respectively, and both models show that the average tensile length tends to a higher value at a rotational speed of 49 r/min, which indicates that higher rotational speeds have vigorous particle movement, which is conducive to the uniform distribution of particles.
Comprehensive analysis of the distribution index and the average stretch length after filtration of the two models leads to the conclusion that the distribution effect is best when the rotational speed is set at 49 r/min. This indicates that increasing the rotational speed is favorable to improving the mixing effect of rubber. This indicates that the 2D and 3D models are consistent in exploring the effect of rubber mixing speed on the mixing effect. Scholar Salahudeen’s study [28] has pointed out that when the speed of rubber mixing machine increases, the tensile length of rubber also increases, and this conclusion is consistent with the conclusion obtained from the study in this paper. The effectiveness of the numerical simulation method is proven.
5.5. Computational Efficiency
The numerical simulation was performed on a Dell personal workstation equipped with two Intel Xeon Platinum 8,269 CY CPUs and 128 GB of RAM. Figure 23 shows a comparison of the computational efficiency of the 2D model and the 3D model, where the horizontal coordinate represents the mixing time of the numerical simulation and the vertical coordinate represents the actual time consumed by the time computer to complete the corresponding numerical simulation. Comparing the computational efficiency of the 2D and 3D models, the computational efficiency of the 2D model is approximately five times that of the 3D model.
6. Conclusions
The numerical simulation method was used to setup three different rotational speed working conditions. By comparing and analyzing the results obtained from 2D and 3D model calculations through evaluation indices such as mixing index, cumulative probability distribution of maximum shear stress, distribution index, and mean length of stretch, the following conclusions can be drawn:(i)The mixing simulation model was established, and the effectiveness of the simulation method was verified by the evaluation and analysis of the mixing results at different speeds and the comparison with the results of open literature.(ii)The speed has a great influence on the mixing. When the speed is set to 49 r/min, it shows the best dispersion and mixing effects and can achieve the effect of uniform distribution faster, showing the best mixing effect.(iii)The changes and trends of the evaluation indicators of the 2D and 3D models are consistent, and the conclusions drawn by the two models are the same. The calculation efficiency of the 2D model is much higher than that of the 3D model.(iv)Within the working speed range of the rubber mixer, increasing the speed can improve the mixing effect of rubber.
It should be noted here that the conditions of this study are isothermal conditions, and the energy equation is not considered. The above conclusions are applicable to isothermal simulation conditions. The applicability of the above conclusions to nonisothermal conditions still needs to be studied. This is the limitation of this study. The advantages and innovations of this research are that the consistency of 2D and 3D models is proved through detailed research, which provides a basis for future scholars to choose 2D and 3D models. However, due to the limitation of the test conditions, the test of the effect of rotational speed on the rubber mixing effect was not conducted. This is the drawback of this study. However, in order to demonstrate the accuracy of the results, it is proved that the conclusions of this research are consistent with the conclusions of existing public researches by consulting the public research of relevant scholars, which has been elaborated in the previous article. On the other hand, it proves the accuracy of the conclusion.
Data Availability
All data, models, and algorithm structures in this paper are available.
Conflicts of Interest
The authors have no conflicts of interest to declare that are relevant to the content of this article.
Acknowledgments
The authors are pleased to acknowledge funding by the National Natural Science Foundation of China (Nos. 52072156, 52272366, 51605198) and the Postdoctoral Foundation of China (2020M682269).