Simulation Studies on the Effect of Material Characteristics and Runners Layout Geometry on the Filling Imbalance in Geometrically Balanced Injection Molds

Simulation studies were performed on filling imbalance in geometrically balanced injection molds. A special simulation procedure was applied to simulate properly the phenomenon, including inertia effects and 3D tetrahedron meshing as well as meshing of the nozzle. The phenomenon was investigated by simulation using several different runner systems at various thermo-rheological material parameters and process operating conditions. It has been observed that the Cross-WLF parameters, index flow, critical shear stress (relaxation time), and zero viscosity, as well as thermal diffusivity and heat transfer coefficient strongly affect the filling imbalance. The effect is substantially dependent on the runners’ layout geometry, as well as on the operating conditions, flow rate, and shear rate. The standard layout geometry and the corrected layout with circled element induce positive imbalance which means that inner cavities fills out faster, and it is opposite for the corrected layouts with one/two overturn elements which cause negative imbalance. Generally, for the standard layout geometry and the corrected layout with circled element, an effect of the zero shear rate viscosity η0 is positive (imbalance increases with an increase of viscosity), and an effect of the power law index n and the relaxation time λ is negative (imbalance decreases with an increase of index n and relaxation time λ). An effect of the thermal diffusivity α and heat transfer coefficient h is negative while an effect of the shear rate is positive. For the corrected layouts with one/two overturn elements, the results of simulations indicate opposite relationships. A novel optimization approach solving the filling imbalance problem and a novel concept of global modeling of injection molding process are also discussed.


Introduction
The phenomenon of filling imbalance in multi-cavity injection molds equipped with geometrically balanced runner systems is a serious problem in the injection molding process. This phenomenon can be observed in H-type runner layouts, e.g., with eight cavities (Figure 1). Additionally, it is generally established that the filling imbalance is caused by the shear gradients in the flow through the runners that lead to non-linear and non-symmetrical distributions of temperature and viscosity. This is strongly influenced and complicated by the geometry of the runners, the thermo-rheological material characteristics, and the process parameters of injection molding [1][2][3][4][5][6]. The filling imbalance has been examined extensively by scientists and engineers, e.g., . Up to now, there have been no universal and commonly accepted solutions of this problem. A detailed discussion of the literature review has been recently presented by the authors of this paper [1].
The filling imbalance has been also studied in the polymer industry. A good example of injection molding process very sensitive to filling imbalance is the manufacturing of lenses [26]. An importance of the effect of filling imbalance on the mold cost and energy consumption has been presented in [27].
The fundamental research on the problem had been carried out by Beaumont group who developed melt rotation technology to minimize the filling imbalance [2][3][4][5]. Beaumont proposed several runner layouts, however, these concepts were not subjected to systematic simulations or experimentations.
Recently, extensive experimental studies on the filling imbalance have been performed by the authors [1]. Balancing the polymer melt flow between cavities has been investigated at various operating conditions using different runner systems. The experiments indicate that the injection rate, the mold and melt temperatures substantially affect the filling imbalance. Additionally, the filling imbalance is strongly dependent on the runners' layout geometry. It is worth noticing that filling imbalance has never been eliminated completely.
The flows of molten polymers in injection molds are unsteady, non-Newtonian and nonisothermal flows occurring at very high rates of deformation. The primary cause of the filling imbalance are nonlinear distributions of the polymer melt velocity and shear rate, as well as the melt temperature which determine the melt viscosity distribution. Such complex flows can only be modeled numerically, e.g., using FEM (finite element method) computations.
First, numerical simulations of the filling imbalance [28][29][30][31] did not allow to predict the phenomenon properly. The use of 2D or 2.5D approaches did not allow to solve the problem. Thus, new concepts were necessary which included inertia effects and 3D, non-isothermal, and non-Newtonian flow [32,33]. Several simulation tests have confirmed the effectiveness of these methods [1,[34][35][36][37].
When modeling the flow of molten polymers in injection molds, the initial flow parameters, e.g., melt temperature and pressure, are basically unknown and usually assumed without much justification. Therefore, a comprehensive approach to modeling the injection molding process is proposed which includes the polymer melt flow both in the plasticating unit as well as in the mold.
Modeling of the polymer melt flow in the plasticating unit of the injection molding machine is based on the modeling of the extrusion process. The first fundamental research in this area has been carried out by Tadmor et al. [38,39] for single screw extrusion. These studies have been later extended to twin screw extrusion, both co-rotating [40,41] and counter-rotating [42][43][44], and starve-fed single screw extrusion [45,46], and finally have been successfully used for modeling of the injection molding process [47][48][49][50]. The modeling research on the screw processing has been summarized by Altinkaynak et al. [51] and Wilczyński et al. [52]. The mechanism of polymer melting in the single screw extruder has been depicted in Figure 2, both for flood fed extrusion (classical Tadmor mechanism) and for starve fed extrusion (two-stage mechanism).
The authors have carried out an experimental (unpublished) research on the melting mechanism of polymers in the injection molding machine (Figure 3). It is interesting to note that the melting in The filling imbalance has been examined extensively by scientists and engineers, e.g., . Up to now, there have been no universal and commonly accepted solutions of this problem. A detailed discussion of the literature review has been recently presented by the authors of this paper [1].
The filling imbalance has been also studied in the polymer industry. A good example of injection molding process very sensitive to filling imbalance is the manufacturing of lenses [26]. An importance of the effect of filling imbalance on the mold cost and energy consumption has been presented in [27].
The fundamental research on the problem had been carried out by Beaumont group who developed melt rotation technology to minimize the filling imbalance [2][3][4][5]. Beaumont proposed several runner layouts, however, these concepts were not subjected to systematic simulations or experimentations.
Recently, extensive experimental studies on the filling imbalance have been performed by the authors [1]. Balancing the polymer melt flow between cavities has been investigated at various operating conditions using different runner systems. The experiments indicate that the injection rate, the mold and melt temperatures substantially affect the filling imbalance. Additionally, the filling imbalance is strongly dependent on the runners' layout geometry. It is worth noticing that filling imbalance has never been eliminated completely.
The flows of molten polymers in injection molds are unsteady, non-Newtonian and non-isothermal flows occurring at very high rates of deformation. The primary cause of the filling imbalance are nonlinear distributions of the polymer melt velocity and shear rate, as well as the melt temperature which determine the melt viscosity distribution. Such complex flows can only be modeled numerically, e.g., using FEM (finite element method) computations.
First, numerical simulations of the filling imbalance [28][29][30][31] did not allow to predict the phenomenon properly. The use of 2D or 2.5D approaches did not allow to solve the problem. Thus, new concepts were necessary which included inertia effects and 3D, non-isothermal, and non-Newtonian flow [32,33]. Several simulation tests have confirmed the effectiveness of these methods [1,[34][35][36][37].
When modeling the flow of molten polymers in injection molds, the initial flow parameters, e.g., melt temperature and pressure, are basically unknown and usually assumed without much justification. Therefore, a comprehensive approach to modeling the injection molding process is proposed which includes the polymer melt flow both in the plasticating unit as well as in the mold.
Modeling of the polymer melt flow in the plasticating unit of the injection molding machine is based on the modeling of the extrusion process. The first fundamental research in this area has been carried out by Tadmor et al. [38,39] for single screw extrusion. These studies have been later extended to twin screw extrusion, both co-rotating [40,41] and counter-rotating [42][43][44], and starve-fed single screw extrusion [45,46], and finally have been successfully used for modeling of the injection molding process [47][48][49][50]. The modeling research on the screw processing has been summarized by Altinkaynak et al. [51] and Wilczyński et al. [52]. The mechanism of polymer melting in the single screw extruder has been depicted in Figure 2, both for flood fed extrusion (classical Tadmor mechanism) and for starve fed extrusion (two-stage mechanism).
The authors have carried out an experimental (unpublished) research on the melting mechanism of polymers in the injection molding machine (Figure 3). It is interesting to note that the melting in the injection molding machine occurs to some extent according to the Tadmor mechanism, however, with clearly visible starvation.
The existing models of injection molding process (plasticating unit) [47][48][49][50] differ from the extrusion models in that they involve the static and dynamic phases of melting (stationary and rotating screw) with an axial screw movement. However, they assume the screw channel is fully filled with a material as in the flood fed extrusion (Figure 2a). According to our observation ( Figure 3) it is not true, and starvation is here clearly seen as in the starve fed extrusion (Figure 2b). Thus, it may be supposed that the solid conveying section and melting section should be modelled similarly to the starve-fed extrusion [45,46], of course accepting the other assumptions valid for the injection molding process.
Polymers 2019, 11,639 3 of 20 the injection molding machine occurs to some extent according to the Tadmor mechanism, however, with clearly visible starvation. The existing models of injection molding process (plasticating unit) [47][48][49][50] differ from the extrusion models in that they involve the static and dynamic phases of melting (stationary and rotating screw) with an axial screw movement. However, they assume the screw channel is fully filled with a material as in the flood fed extrusion (Figure 2a). According to our observation ( Figure 3) it is not true, and starvation is here clearly seen as in the starve fed extrusion (Figure 2b). Thus, it may be supposed that the solid conveying section and melting section should be modelled similarly to the starve-fed extrusion [45,46], of course accepting the other assumptions valid for the injection molding process.  Recently, novel concepts of modeling of polymer melting have been presented [51,53]. These involve the modeling of a two-phase flow (solid/melt) as a single-phase flow, and the resulting temperature distribution determines the boundary between melt and solid. This concept is independent on the screw geometry, and does not require any knowledge of the melting mechanism. More advanced concepts include DEM modeling (discrete element method) for solid conveying and FEM modeling (finite element method) for melting and melt flow. It is also worth highlighting that the full 3D simulations are well developed for die extrusion flows including viscoelastic models [54].
So far, the filling imbalance has been mainly studied experimentally and in a much less extent by simulation (for the standard systems and for the basic melt rotation systems, only). In this paper, an extensive simulation study on the filling imbalance has been performed for several different runner systems, at various thermo-rheological material parameters and process operating parameters, to explain the phenomenon. The Autodesk Simulation Moldflow Insight 2014 software has been used for simulations [26]. According to our knowledge it is the most comprehensive simulation and modeling research on the filling imbalance with particular emphasis on material and geometrical effects.
The novel approaches solving the imbalance problem have been also presented including an optimization concept with the use of advanced process modeling, as well as the concept of global modeling of the injection molding process. the injection molding machine occurs to some extent according to the Tadmor mechanism, however, with clearly visible starvation. The existing models of injection molding process (plasticating unit) [47][48][49][50] differ from the extrusion models in that they involve the static and dynamic phases of melting (stationary and rotating screw) with an axial screw movement. However, they assume the screw channel is fully filled with a material as in the flood fed extrusion (Figure 2a). According to our observation ( Figure 3) it is not true, and starvation is here clearly seen as in the starve fed extrusion (Figure 2b). Thus, it may be supposed that the solid conveying section and melting section should be modelled similarly to the starve-fed extrusion [45,46], of course accepting the other assumptions valid for the injection molding process.  Recently, novel concepts of modeling of polymer melting have been presented [51,53]. These involve the modeling of a two-phase flow (solid/melt) as a single-phase flow, and the resulting temperature distribution determines the boundary between melt and solid. This concept is independent on the screw geometry, and does not require any knowledge of the melting mechanism. More advanced concepts include DEM modeling (discrete element method) for solid conveying and FEM modeling (finite element method) for melting and melt flow. It is also worth highlighting that the full 3D simulations are well developed for die extrusion flows including viscoelastic models [54].

Process Simulation
So far, the filling imbalance has been mainly studied experimentally and in a much less extent by simulation (for the standard systems and for the basic melt rotation systems, only). In this paper, an extensive simulation study on the filling imbalance has been performed for several different runner systems, at various thermo-rheological material parameters and process operating parameters, to explain the phenomenon. The Autodesk Simulation Moldflow Insight 2014 software has been used for simulations [26]. According to our knowledge it is the most comprehensive simulation and modeling research on the filling imbalance with particular emphasis on material and geometrical effects.
The novel approaches solving the imbalance problem have been also presented including an optimization concept with the use of advanced process modeling, as well as the concept of global modeling of the injection molding process. Recently, novel concepts of modeling of polymer melting have been presented [51,53]. These involve the modeling of a two-phase flow (solid/melt) as a single-phase flow, and the resulting temperature distribution determines the boundary between melt and solid. This concept is independent on the screw geometry, and does not require any knowledge of the melting mechanism. More advanced concepts include DEM modeling (discrete element method) for solid conveying and FEM modeling (finite element method) for melting and melt flow. It is also worth highlighting that the full 3D simulations are well developed for die extrusion flows including viscoelastic models [54].

Process Simulation
So far, the filling imbalance has been mainly studied experimentally and in a much less extent by simulation (for the standard systems and for the basic melt rotation systems, only). In this paper, an extensive simulation study on the filling imbalance has been performed for several different runner systems, at various thermo-rheological material parameters and process operating parameters, to explain the phenomenon. The Autodesk Simulation Moldflow Insight 2014 software has been used for simulations [26]. According to our knowledge it is the most comprehensive simulation and modeling research on the filling imbalance with particular emphasis on material and geometrical effects.
The novel approaches solving the imbalance problem have been also presented including an optimization concept with the use of advanced process modeling, as well as the concept of global modeling of the injection molding process.

Simulation Program
Process simulations have been carried out to study an effect of the material characteristics and runners' layout geometry on the filling imbalance in geometrically balanced injection molds.
Simulations have been performed for an eight-cavity injection mold of "H-type" runners' layout equipped with inserts of different geometry. A standard geometry GS has been used, as well as three overturn geometries: G1 with one correction element, G2 with two correction elements, and G3 with circled element, These are depicted in Figure 4.
We have used polybutylene terephthalate PBT Valox 337 in the study (manufactured by SABIC, Riyadh, Saudi Arabia) as the reference polymer to compare the research results with recently performed experimentations [1] and other simulations [31,32]. The material characteristics has been collected in Table 1 [28].
The viscosity of the material was modeled using the non-Newtonian Cross-WLF model which describes the shear thinning behaviour and the temeperature effect, and it may be presented in the following form: where is the melt viscosity, Pa•s; is the zero viscosity, Pa•s; is the shear rate, s −1 ; * is the critical shear stress, Pa; / * is the relaxation time λ, s; n is the power law index, dimensionless; T is the temperature, K; and D1, A1, A2, and T * are the parameters of the Cross-WLF equation.   An effect of the thermo-rheological material characteristics has been studied under different operating conditions, at various injection rates (flow rates). The processing parameters are presented in Table 1.

Material
We have used polybutylene terephthalate PBT Valox 337 in the study (manufactured by SABIC, Riyadh, Saudi Arabia) as the reference polymer to compare the research results with recently performed experimentations [1] and other simulations [31,32]. The material characteristics has been collected in Table 1 [28].
The viscosity of the material was modeled using the non-Newtonian Cross-WLF model which describes the shear thinning behaviour and the temeperature effect, and it may be presented in the following form: where η is the melt viscosity, Pa·s; η 0 is the zero viscosity, Pa·s; . γ is the shear rate, s −1 ; τ * is the critical shear stress, Pa; η 0 /τ * is the relaxation time λ, s; n is the power law index, dimensionless; T is the temperature, K; and D 1 , A 1 , A 2 , and T * are the parameters of the Cross-WLF equation.
An effect of the Cross-WLF parameters on the shear viscosity curves is depicted in Figures 5-7 in relation to the parameters of the reference material ( Table 1).
The two limiting values of the Cross-WLF model parameters have been selected as the limiting values of these parameters of the materials available in the Autodesk Moldflow database. The different melt temperature was set to change the zero shear viscosity, and the relaxation time λ was set by changing τ*. An effect of the Cross-WLF parameters on the shear viscosity curves is depicted in Figures 5-7 in relation to the parameters of the reference material ( Table 1).
The two limiting values of the Cross-WLF model parameters have been selected as the limiting values of these parameters of the materials available in the Autodesk Moldflow database. The different melt temperature was set to change the zero shear viscosity, and the relaxation time λ was set by changing τ*.
A sharp drop of the shear viscosity for low values of the power law index n (shear thinning) at higher shear rates (typical for injection molding) is clearly visible ( Figure 5). An increase in the relaxation time λ also leads to the significant lowering of the viscosity at higher shear rates . An effect of the zero shear viscosity on the viscosity curve is greater for lower shear rates , however, also significant in the injection molding shear rate range. In summary, it can be stated that rheological parameters of the material have a significant impact on the shear viscosity in the typical for injection molding range of the shear rate, that is = (10 3 -10 4 , 10 5 ) s −1 . This affects the filling imbalance phenomenon.

Procedure of Estimation of Filling Imbalance
The mass filling imbalance coefficient Im [2] has been used to estimate the filling imbalance which is defined by: where Im is the mass filling imbalance coefficient, m1 is the mass of the polymer in the inner cavities, and m2 is the mass of the polymer in the outer cavities ( Figure 8).
The mass filling imbalance coefficient is positive when the inner cavities (1) are filled out faster (Im > 0), and it is negative when the outer cavities (2) are filled out faster (Im < 0).

Procedure of Estimation of Filling Imbalance
The mass filling imbalance coefficient Im [2] has been used to estimate the filling imbalance which is defined by: where Im is the mass filling imbalance coefficient, m1 is the mass of the polymer in the inner cavities, and m2 is the mass of the polymer in the outer cavities ( Figure 8).
The mass filling imbalance coefficient is positive when the inner cavities (1) are filled out faster (Im > 0), and it is negative when the outer cavities (2) are filled out faster (Im < 0).  A sharp drop of the shear viscosity η for low values of the power law index n (shear thinning) at higher shear rates . γ (typical for injection molding) is clearly visible ( Figure 5). An increase in the relaxation time λ also leads to the significant lowering of the viscosity η at higher shear rates .

Simulation Technique
γ. An effect of the zero shear viscosity η 0 on the viscosity curve is greater for lower shear rates . γ, however, also significant in the injection molding shear rate range. In summary, it can be stated that rheological parameters of the material have a significant impact on the shear viscosity η in the typical for injection molding range of the shear rate, that is . γ = (10 3 -10 4 , 10 5 ) s −1 . This affects the filling imbalance phenomenon.

Procedure of Estimation of Filling Imbalance
The mass filling imbalance coefficient I m [2] has been used to estimate the filling imbalance which is defined by: where I m is the mass filling imbalance coefficient, m 1 is the mass of the polymer in the inner cavities, and m 2 is the mass of the polymer in the outer cavities ( Figure 8). The mass filling imbalance coefficient is positive when the inner cavities (1) are filled out faster (I m > 0), and it is negative when the outer cavities (2) are filled out faster (I m < 0).

1
100% (3) where Im is the mass filling imbalance coefficient, m1 is the mass of the polymer in the inner cavities, and m2 is the mass of the polymer in the outer cavities ( Figure 8). The mass filling imbalance coefficient is positive when the inner cavities (1) are filled out faster (Im > 0), and it is negative when the outer cavities (2) are filled out faster (Im < 0).

Simulation Technique
Simulations have been performed using the Moldflow (San Rafael, CA, USA) Insight 2014 program [28]. A special simulation technique has been applied to study the filling imbalance phenomenon. It embraces the geometrical model design of the flow space, meshing technique and definition of the solver parameters: -the geometrical model of the flow space should cover all the flow and the flow environment components, i.e., the cavities, the runners, the moldbase, and the cooling channels, as well as the nozzle where the filling imbalance starts to develop (Figure 9a), and the proper design of cooling channels ensures an even temperature distribution of the mold (Figure 9b).
the meshing should be performed using the 3D tetrahedron mesh elements to properly model the asymmetric distributions of the flow parameters, e.g., shear rate and temperature, and the maksimum mesh size should not exceed the f value (f = D/2N, D-characteristic dimension of the channel, cavity or runners, N-number of layers), and the number of layers across the part thickness should not be lower than 12 (the only fine solid mesh is able to catch the shear rate distribution properly across the runner), -when modeling and using the equation of motion, the inertia components should be taken into account which improves the prediction of the polymer melt front. Simulations have been performed using the Moldflow (San Rafael, CA, USA) Insight 2014 program [28]. A special simulation technique has been applied to study the filling imbalance phenomenon. It embraces the geometrical model design of the flow space, meshing technique and definition of the solver parameters: -the geometrical model of the flow space should cover all the flow and the flow environment components, i.e., the cavities, the runners, the moldbase, and the cooling channels, as well as the nozzle where the filling imbalance starts to develop (Figure 9a), and the proper design of cooling channels ensures an even temperature distribution of the mold (Figure 9b).
the meshing should be performed using the 3D tetrahedron mesh elements to properly model the asymmetric distributions of the flow parameters, e.g., shear rate and temperature, and the maksimum mesh size should not exceed the f value (f = D/2N, D-characteristic dimension of the channel, cavity or runners, N-number of layers), and the number of layers across the part thickness should not be lower than 12 (the only fine solid mesh is able to catch the shear rate distribution properly across the runner), -when modeling and using the equation of motion, the inertia components should be taken into account which improves the prediction of the polymer melt front. It is important to know, that the modeling of filling imbalance is time consuming and usually takes several hours (6)(7)(8) for about 1,500,000 tetrahedrons. The simulations were performed using an HP Z640 workstation; Intel ® Xeon ® E5-2667 v4 (3.2 GHz, 25 MB cache, 8 cores, Intel ® vPro™) equipped with 128 GB DDR4-2400 ECC registered RAM.
It is also important to know that the conventional 1D beam elements are not suitable for the runner modeling (Figure 10a). These cannot catch the non-symmetrical temperature distribution in the runner cross section. Combination of the conventional 1D beam elements with 3D elements is also not recommended. At the transition cross-section between these two meshes the shear rate parameters are not transferred properly (Figure 10b). Only the 3D mesh is able to simulate the shear induced filling imbalance properly (Figure 10c). It is important to know, that the modeling of filling imbalance is time consuming and usually takes several hours (6)(7)(8) for about 1,500,000 tetrahedrons. The simulations were performed using an HP Z640 workstation; Intel ® Xeon ® E5-2667 v4 (3.2 GHz, 25 MB cache, 8 cores, Intel ® vPro™) equipped with 128 GB DDR4-2400 ECC registered RAM.
It is also important to know that the conventional 1D beam elements are not suitable for the runner modeling (Figure 10a). These cannot catch the non-symmetrical temperature distribution in the runner cross section. Combination of the conventional 1D beam elements with 3D elements is also not recommended. At the transition cross-section between these two meshes the shear rate parameters are not transferred properly (Figure 10b). Only the 3D mesh is able to simulate the shear induced filling imbalance properly (Figure 10c).

Results
Simulations have been performed for all the runner layouts (Figure 4), i.e., for a standard geometry and for three overturn geometries, at various injection rates Vinj, for different rheological and thermal parameters of the material. An effect of the parameters of the Cross-WLF model, i.e., the zero shear rate viscosity η0, the power law index n and the relaxation time λ, has been studied as well as an effect of the thermal properties of the material, i.e., the thermal diffusivity α and the heat transfer coefficient h has been simulated. Obviously, while considering the effect of the selected parameter, the other were kept constant.
Mass filling imbalance coefficient Im has been plotted against injection rate Vinj, and distributions of the polymer melt velocity have been depicted for all the simulation cases.
The surface areas of the filled parts of cavities have been measured at the stage of 80% of mold filling, and these areas have been inserted in places of m1 and m2 (Equation 3). It has been assumed the density is constant.
First, an effect of the rheological parameters has been studied, then an effect of the thermal parameters has been discussed and, finally, an effect of the runner geometry for selected rheological and thermal parameters has been considered.

Effect of Rheological Parameters
An effect of the rheological parameters of the material on the filling imbalance has been depicted in Figures 11-16. Here, as an example, the simulation results for the standard geometry of the runners' layout have been presented for the lowest and the highest values of the parameters studied ( Figures 5-7).

The Zero Shear Viscosity
An effect of the zero shear rate viscosity η0 on the filling imbalance for the standard geometry of

Results
Simulations have been performed for all the runner layouts (Figure 4), i.e., for a standard geometry and for three overturn geometries, at various injection rates V inj , for different rheological and thermal parameters of the material. An effect of the parameters of the Cross-WLF model, i.e., the zero shear rate viscosity η 0 , the power law index n and the relaxation time λ, has been studied as well as an effect of the thermal properties of the material, i.e., the thermal diffusivity α and the heat transfer coefficient h has been simulated. Obviously, while considering the effect of the selected parameter, the other were kept constant.
Mass filling imbalance coefficient I m has been plotted against injection rate V inj , and distributions of the polymer melt velocity have been depicted for all the simulation cases.
The surface areas of the filled parts of cavities have been measured at the stage of 80% of mold filling, and these areas have been inserted in places of m 1 and m 2 (Equation (3)). It has been assumed the density is constant.
First, an effect of the rheological parameters has been studied, then an effect of the thermal parameters has been discussed and, finally, an effect of the runner geometry for selected rheological and thermal parameters has been considered.

Effect of Rheological Parameters
An effect of the rheological parameters of the material on the filling imbalance has been depicted in Figures 11-16. Here, as an example, the simulation results for the standard geometry of the runners' layout have been presented for the lowest and the highest values of the parameters studied ( Figures 5-7).

. The Zero Shear Viscosity
An effect of the zero shear rate viscosity η 0 on the filling imbalance for the standard geometry of the runners layout is generally positive which means that the inner cavities fills out faster, and this significantly increases with an increase of the injection rate ( Figure 11). The imbalance increases when the viscosity increases. Higher viscosity means the higher energy dissipation, so it means the flow with higher temperature gradients which promotes the filling imbalance. Velocity distributions and the flow front advancement (Figure 12) clearly confirm this observation, and the higher velocity gradients are observed when the viscosity increases.

The Relaxation Time
An effect of the relaxation time λ on the filling imbalance is positive for low λ values, and then the inner cavities fills out faster, and this increases significantly with an increase of the injection rate ( Figure 13). However, the imbalance decreases almost to zero when the relaxation time increases. This results from a substantial viscosity decrease (Figure 6), i.e., the lower energy dissipation, and the flow with lower temperature gradients which does not promote the filling imbalance. Velocity distributions and the flow front advancement (Figure 14) confirm this observation, and the higher velocity gradients are observed for the lower relaxation time.

The Relaxation Time
An effect of the relaxation time λ on the filling imbalance is positive for low λ values, and then the inner cavities fills out faster, and this increases significantly with an increase of the injection rate ( Figure 13). However, the imbalance decreases almost to zero when the relaxation time increases. This results from a substantial viscosity decrease (Figure 6), i.e., the lower energy dissipation, and the flow with lower temperature gradients which does not promote the filling imbalance. Velocity distributions and the flow front advancement (Figure 14) confirm this observation, and the higher velocity gradients are observed for the lower relaxation time.

The Relaxation Time
An effect of the relaxation time λ on the filling imbalance is positive for low λ values, and then the inner cavities fills out faster, and this increases significantly with an increase of the injection rate ( Figure 13). However, the imbalance decreases almost to zero when the relaxation time increases. This results from a substantial viscosity decrease (Figure 6), i.e., the lower energy dissipation, and the flow with lower temperature gradients which does not promote the filling imbalance. Velocity distributions and the flow front advancement (Figure 14) confirm this observation, and the higher velocity gradients are observed for the lower relaxation time.

The Power Law Index
An effect of the power law index n is also positive for low n values which means the inner cavities fills out faster, and this increases substantially with an increase of the injection rate ( Figure 15). However, the imbalance decreases basically to zero when the power low index increases. This is a clear proof that the imbalance is greater when the fluid is more non-Newtonian. Velocity distributions and the flow front advancement (Figure 16) confirm this observation, and the higher velocity gradients are observed for the lower power law index.

The Relaxation Time
An effect of the relaxation time λ on the filling imbalance is positive for low λ values, and then the inner cavities fills out faster, and this increases significantly with an increase of the injection rate ( Figure 13). However, the imbalance decreases almost to zero when the relaxation time increases. This results from a substantial viscosity decrease (Figure 6), i.e., the lower energy dissipation, and the flow with lower temperature gradients which does not promote the filling imbalance. Velocity distributions and the flow front advancement (Figure 14) confirm this observation, and the higher velocity gradients are observed for the lower relaxation time.

The Power Law Index
An effect of the power law index n is also positive for low n values which means the inner cavities fills out faster, and this increases substantially with an increase of the injection rate ( Figure 15). However, the imbalance decreases basically to zero when the power low index increases. This is a clear proof that the imbalance is greater when the fluid is more non-Newtonian. Velocity distributions and the flow front advancement (Figure 16) confirm this observation, and the higher velocity gradients are observed for the lower power law index.

The Power Law Index
An effect of the power law index n is also positive for low n values which means the inner cavities fills out faster, and this increases substantially with an increase of the injection rate ( Figure 15). However, the imbalance decreases basically to zero when the power low index increases. This is a clear proof that the imbalance is greater when the fluid is more non-Newtonian. Velocity distributions and the flow front advancement (Figure 16) confirm this observation, and the higher velocity gradients are observed for the lower power law index.

Effect of Thermal Parameters
An effect of the thermal parameters of the material on the filling imbalance has been depicted in Figures 17-20. Here, as an example, the simulation results for the standard geometry of the runners' layout have been presented.

The Thermal Diffusivity
The thermal diffusivity α (the middle value α = 0.112 mm 2 /s) which is a combination of the thermal conductivity, density, and specific heat capacity, was set by changing the thermal conductivity.
An effect of the thermal diffusivity α on the filling imbalance is positive for low and high α values which means that inner cavities fills out faster, This positive effect significantly increases with an increase of the injection rate ( Figure 17). However, the imbalance decreases when the thermal diffusivity increases. The results of the simulation show that a more intense heat exchange, i.e., higher α promotes balancing the flow. This is clearly confirmed by velocity distributions and the flow front advancement (Figure 18). clear proof that the imbalance is greater when the fluid is more non-Newtonian. Velocity distributions and the flow front advancement (Figure 16) confirm this observation, and the higher velocity gradients are observed for the lower power law index.

Effect of Thermal Parameters
An effect of the thermal parameters of the material on the filling imbalance has been depicted in Figures 17-20. Here, as an example, the simulation results for the standard geometry of the runners' layout have been presented.

The Thermal Diffusivity
The thermal diffusivity α (the middle value α = 0.112 mm 2 /s) which is a combination of the thermal conductivity, density, and specific heat capacity, was set by changing the thermal conductivity.
An effect of the thermal diffusivity α on the filling imbalance is positive for low and high α values which means that inner cavities fills out faster, This positive effect significantly increases with an increase of the injection rate ( Figure 17). However, the imbalance decreases when the thermal diffusivity increases. The results of the simulation show that a more intense heat exchange, i.e., higher α promotes balancing the flow. This is clearly confirmed by velocity distributions and the flow front advancement (Figure 18).

Effect of Thermal Parameters
An effect of the thermal parameters of the material on the filling imbalance has been depicted in Figures 17-20. Here, as an example, the simulation results for the standard geometry of the runners' layout have been presented.

The Thermal Diffusivity
The thermal diffusivity α (the middle value α = 0.112 mm 2 /s) which is a combination of the thermal conductivity, density, and specific heat capacity, was set by changing the thermal conductivity.
An effect of the thermal diffusivity α on the filling imbalance is positive for low and high α values which means that inner cavities fills out faster, This positive effect significantly increases with an increase of the injection rate ( Figure 17). However, the imbalance decreases when the thermal diffusivity increases. The results of the simulation show that a more intense heat exchange, i.e., higher α promotes balancing the flow. This is clearly confirmed by velocity distributions and the flow front advancement (Figure 18).

The Heat Transfer Coefficient
The heat transfer coefficient describes the heat transfer between the mold and the polymer melt. The middle value h = 5000 W/(m 2 K) has been assumed by default (Autodesk Moldflow).
An effect of the heat transfer coefficient h on the filling imbalance is also positive which means

The Heat Transfer Coefficient
The heat transfer coefficient describes the heat transfer between the mold and the polymer melt. The middle value h = 5000 W/(m 2 K) has been assumed by default (Autodesk Moldflow).
An effect of the heat transfer coefficient h on the filling imbalance is also positive which means that inner cavities fills out faster. However, the relationship between the imbalance, the heat transfer coefficient, and the flow rate is not very clear (Figure 19). For low h values, the imbalance is higher for the low flow rate, while for high h values, the imbalance is higher for the high flow rate. This is confirmed by Figure 20. Additionally, it is worth noting that the selection of an appropriate value of the h coefficient is important to make the phenomenon visible by simulation. Unfortunatelly, there is generally a lack of good material data in this area.
Polymers 2019, 11,639 12 of 20 coefficient, and the flow rate is not very clear ( Figure 19). For low h values, the imbalance is higher for the low flow rate, while for high h values, the imbalance is higher for the high flow rate. This is confirmed by Figure 20. Additionally, it is worth noting that the selection of an appropriate value of the h coefficient is important to make the phenomenon visible by simulation. Unfortunatelly, there is generally a lack of good material data in this area.

Geometry Effects
An effect of the runner geometry for the selected rheological and thermal material parameters has been studied. The relaxation time λ and the thermal diffusivity α have been chosen for this. Here, as an example, the simulation results for the standard geometry GS of the runners layout, as well for the one overturn geometry G1 have been presented for the high flow rate (500 cm 3 /s) (Figures 21-24).

The Relaxation Time
For the standard geometry GS (Figure 4a, Figure 21), the filling imbalance is positive (Im = 20%) for the low λ value (λ = 0.001 s) which means that inner cavities fills out faster, However, the imbalance decreases almost to zero (Im = 1%), when the relaxation time increases (λ = 0.1 s). This is consistent with Figure 13.
For the one overturn geometry G1 (Figure 4b, Figure 22), the filling imbalance is negative (Im = −35%) for the low λ value (λ = 0.001 s) which means that outer cavities fills out faster, but the imbalance decreases basically to zero (Im = 1%), when the relaxation time increases (λ = 0.1 s).
Thus, in both cases, the geometry GS and G1, the imbalance is very high for the low relaxation coefficient, and the flow rate is not very clear ( Figure 19). For low h values, the imbalance is higher for the low flow rate, while for high h values, the imbalance is higher for the high flow rate. This is confirmed by Figure 20. Additionally, it is worth noting that the selection of an appropriate value of the h coefficient is important to make the phenomenon visible by simulation. Unfortunatelly, there is generally a lack of good material data in this area.

Geometry Effects
An effect of the runner geometry for the selected rheological and thermal material parameters has been studied. The relaxation time λ and the thermal diffusivity α have been chosen for this. Here, as an example, the simulation results for the standard geometry GS of the runners layout, as well for the one overturn geometry G1 have been presented for the high flow rate (500 cm 3 /s) (Figures 21-24).

The Relaxation Time
For the standard geometry GS (Figure 4a, Figure 21), the filling imbalance is positive (Im = 20%) for the low λ value (λ = 0.001 s) which means that inner cavities fills out faster, However, the imbalance decreases almost to zero (Im = 1%), when the relaxation time increases (λ = 0.1 s). This is consistent with Figure 13.
For the one overturn geometry G1 (Figure 4b, Figure 22), the filling imbalance is negative (Im = −35%) for the low λ value (λ = 0.001 s) which means that outer cavities fills out faster, but the imbalance decreases basically to zero (Im = 1%), when the relaxation time increases (λ = 0.1 s).
Thus, in both cases, the geometry GS and G1, the imbalance is very high for the low relaxation time, and almost disappears for the high relaxation time. This results from a significant viscosity

Geometry Effects
An effect of the runner geometry for the selected rheological and thermal material parameters has been studied. The relaxation time λ and the thermal diffusivity α have been chosen for this. Here, as an example, the simulation results for the standard geometry GS of the runners layout, as well for the one overturn geometry G1 have been presented for the high flow rate (500 cm 3 /s) (Figures 21-24).

The Relaxation Time
For the standard geometry GS (Figures 4a and 21), the filling imbalance is positive (I m = 20%) for the low λ value (λ = 0.001 s) which means that inner cavities fills out faster, However, the imbalance decreases almost to zero (I m = 1%), when the relaxation time increases (λ = 0.1 s). This is consistent with Figure 13. For the one overturn geometry G1 (Figures 4b and 22), the filling imbalance is negative (I m = −35%) for the low λ value (λ = 0.001 s) which means that outer cavities fills out faster, but the imbalance decreases basically to zero (I m = 1%), when the relaxation time increases (λ = 0.1 s).
Thus, in both cases, the geometry GS and G1, the imbalance is very high for the low relaxation time, and almost disappears for the high relaxation time. This results from a significant viscosity decrease with a relaxation time increase (Figure 6), i.e., the lower energy dissipation. This causes the flow with lower temperature gradients which does not promote the filling imbalance.
It is clearly seen that the high imbalance results from the high temperature gradients (Figures 21a  and 22a), and when the flow is balanced, the temperature gradients disappear (Figures 21b and 22b).
It is also clearly seen that for the standard geometry GS the polymer melt stream rotates to the right (cross-section C-C, Figure 21a) which results in a positive imbalance. While, for the one overturn geometry G1 the polymer melt stream rotates to the left (cross-section C-C, Figure 22a) which results in a negative imbalance. It is worth noting that for the high relaxation time (Figures 21a and 22b) rotation is not observed. decrease with a relaxation time increase (Figure 6), i.e., the lower energy dissipation. This causes the flow with lower temperature gradients which does not promote the filling imbalance. It is clearly seen that the high imbalance results from the high temperature gradients ( Figure  21a, Figure 22a), and when the flow is balanced, the temperature gradients disappear (Figure 21b,  Figure 22b).
It is also clearly seen that for the standard geometry GS the polymer melt stream rotates to the right (cross-section C-C, Figure 21a) which results in a positive imbalance. While, for the one overturn geometry G1 the polymer melt stream rotates to the left (cross-section C-C, Figure 22a) which results in a negative imbalance. It is worth noting that for the high relaxation time (Figure 21a, Figure 22b) rotation is not observed.

The Thermal Diffusivity
For the standard geometry GS (Figure 4a, Figure 23), the filling imbalance is positive (Im = 22%) for the low α value (α = 0.005 mm 2 /s) which means that inner cavities fills out faster. Additionally, it is also positive, but lower (Im = 12%), for the high α value (α = 5 mm 2 /s). This is consistent with Figure  17.
In both cases, the geometry GS and G1, the imbalance is higher for the low thermal diffusivity which results from the less intense heat exchange. The more intense heat exchange, i.e., the higher α promotes balancing the flow.

The Thermal Diffusivity
For the standard geometry GS (Figures 4a and 23), the filling imbalance is positive (I m = 22%) for the low α value (α = 0.005 mm 2 /s) which means that inner cavities fills out faster. Additionally, it is also positive, but lower (I m = 12%), for the high α value (α = 5 mm 2 /s). This is consistent with Figure 17.
For the one overturn geometry G1 (Figures 4b and 24), the filling imbalance is negative (I m = −33%) for the low α value (α = 0.005 mm 2 /s) which means that outer cavities fills out faster. Additionally, it is also negative, but lower (I m = 12%), for the high α value (α = 5 mm 2 /s).
In both cases, the geometry GS and G1, the imbalance is higher for the low thermal diffusivity which results from the less intense heat exchange. The more intense heat exchange, i.e., the higher α promotes balancing the flow.
It is clearly seen that the high imbalance results from the high temperature gradients (Figures 23a  and 24a), and the imbalance decreases when lowering the temperature gradients (Figures 23b and 24b).
It is also clearly seen that for the standard geometry GS the polymer melt stream rotates to the right in both cases, the low and high thermal diffusivity (cross-section C-C, Figure 23a,b) which results in a positive imbalance. While, for the one overturn geometry G1 the polymer melt stream rotates to the left (cross-section C-C, Figure 24a,b) which results in a negative imbalance. It is worth noting that rotation is observed in all the cases, which clearly show that imbalance is connected with the temperature gradients and the melt stream rotation. It is clearly seen that the high imbalance results from the high temperature gradients ( Figure  23a, Figure 24a), and the imbalance decreases when lowering the temperature gradients (Figure 23b,  Figure 24b).
It is also clearly seen that for the standard geometry GS the polymer melt stream rotates to the right in both cases, the low and high thermal diffusivity (cross-section C-C, Figure 23a,b) which results in a positive imbalance. While, for the one overturn geometry G1 the polymer melt stream rotates to the left (cross-section C-C, Figure 24a,b) which results in a negative imbalance. It is worth noting that rotation is observed in all the cases, which clearly show that imbalance is connected with the temperature gradients and the melt stream rotation.

Conclusions
It has been observed that the Cross-WLF parameters, index flow, critical shear stress (relaxation time) and zero viscosity, as well as the thermal diffusivity and heat transfer coefficient strongly affect the filling imbalance. The effect is substantially dependent on the runners' layout geometry as well as on the operating conditions, flow rate, and shear rate. Additionally, there are no simple rules for runners' layout design. Generally, the standard runners system GS and the system with circled element G3 induce a positive imbalance which means that inner cavities fills out faster, and it is opposite for the systems with one/two overturn elements G1/G2 which cause negative imbalance. In general, for GS and G3 systems, an effect of the zero shear rate viscosity η0 is positive (imbalance increases with an increase of viscosity), and an effect of the power law index n and the relaxation time λ is negative (imbalance decreases with an increase of index n and relaxation time λ). An effect of the thermal diffusivity α and heat transfer coefficient h is negative while an effect of the shear rate is

Conclusions
It has been observed that the Cross-WLF parameters, index flow, critical shear stress (relaxation time) and zero viscosity, as well as the thermal diffusivity and heat transfer coefficient strongly affect the filling imbalance. The effect is substantially dependent on the runners' layout geometry as well as on the operating conditions, flow rate, and shear rate. Additionally, there are no simple rules for runners' layout design. Generally, the standard runners system GS and the system with circled element G3 induce a positive imbalance which means that inner cavities fills out faster, and it is opposite for the systems with one/two overturn elements G1/G2 which cause negative imbalance. In general, for GS and G3 systems, an effect of the zero shear rate viscosity η 0 is positive (imbalance increases with an increase of viscosity), and an effect of the power law index n and the relaxation time λ is negative (imbalance decreases with an increase of index n and relaxation time λ). An effect of the thermal diffusivity α and heat transfer coefficient h is negative while an effect of the shear rate is positive. For G1 and G2 systems, the results of simulations indicate opposite dependencies. These observations have been summarized in Table 2. Table 2. Simulation summary (green -I m > 0, red -I m < 0).

Effect of Process Parameters and Runners Layout on Filling Imbalance -Summary
It can be concluded, the velocity/shear rate distribution and the heat transfer conditions from the polymer melt to the mold are two main factors affecting the filling imbalance. It was observed that when flow rate increases, the filling imbalance increases since the non-linearity of velocity/shear rate distribution also increases. At low flow rates there is more time for the heat transfer from the polymer melt to the mold, and this may explain the imbalance increase with a decrease of the flow rate.
It can be summarized that filling imbalance problem is still unsolved. There are no universal runner layouts that could be successfully applied for mold design, and the design of runner systems is always dependent on the material characteristics and process parameters. If the design parameters of the runner system are inappropriate, the imbalance still happens. Hence, the assistance of simulation/optimization tools is desired to handle this phenomenon effectively.
The negative effects of mold filling imbalance are widely known in the industry, however, deterministic strategies of reducing this phenomenon are not available. This leads to the conclusion that only the statistical or AI (artificial intelligence) optimization techniques can improve the mold designing. We propose to use the STASA QC [55] system to optimize the filling imbalance applying the brain construction algorithm (BCA) to define a response surface of the problem. The melt injection rate, melt temperature, mold temperature, and the runner geometry will be optimized, and the global objective function will be defined by the filling imbalance coefficient (of the highest weight), injection pressure, molding temperature, and injection time. This concept will be presented in the next paper.
A comprehensive approach to modeling of the injection molding may be also useful for simulation of the flow in injection molds and for the prediction of filling imbalance. A global injection molding model might be considered for simulation of the polymer melt flow in the plasticating unit as well as in the mold. Resulting parameters of the plasticating unit simulations would be input data for the mold flow simulations.