Influences of the mesh in the CAE simulation for plastic injection molding

1Centro de Fabricación Avanzada Aeronáutica – CFAA, Departamento de Ingeniería Mecánica, Universidad del País Vasco/Euskal Herriko Unibertsitatea – UPV/EHU, Bilbao, Vizcay, Espanã 2Grupo de Pesquisa de Manufatura Auxiliada por Computador – GPCAM, Núcleo de Inovação em Moldagem e Manufatura Aditiva – NIMMA, Programa de Pós-graduação em Engenharia Mecânica, Universidade Federal de Santa Catarina – UFSC, Florianópolis, SC, Brasil


Introduction
Plastics have been used for many years to manufacture technical products to supply the automotive, aeronautic, medical and electrometric industries. Injection molding is the main plastic transformation process and several parameters influence the product quality and the cycle time, such as temperature, pressure and cooling time. In many cases, the outcome of the injection molding process is difficult to predict and the definition of the most suitable parameters is obtained empirically. Nowadays, computer-aided engineering (CAE) software can be used to assist the plastics industry, to reduce the costs and the cycle time and improving the product quality [1,2] by assisting in two of the production phases: i) Simulating the plastic product: In CAE, defects in the plastic parts, such as, welding lines and shrinkage, can be simulated [3] . CAE can also be used to detect critical regions with heat accumulation, sink marks, residual stress built-up and product warpage [4,5] . In such cases, the accuracy of the CAE simulation can be easily verified by comparison, accessing the distortion of the manufactured plastic part using a CAD/CAI/ 3D scanner or microscopy for structural analysis and measuring the residual stress. Padilla et al. [6] evaluated the warpage and shrinkage of a product and compared the alterations observed in the experimental injected part with the simulation results. The simulation resulted in some differences in the regions close to edges of the product.
ii) Simulating the molding process: To evaluate the CAE simulation of the plastic injection molding process, an adequate data acquisition system has to be developed and sensors need to be installed inside the mold cavity to obtain the experimental behavior of the plastic during the molding phase. Thus, the accuracy of this type of simulation is not addressed in depth in the current literature In the CAE software, a mesh is used to perform the calculations. The mesh geometry can vary and the mesh can be generated with different densities (number of elements). Theoretically, the higher the mesh density the higher the accuracy of the simulation will be, because there will be more elements to describe the simulated part. However, a longer computational time will be required. Identifying the most appropriate mesh geometry and density to perform the CAE calculation still represents a challenge. The CAE software user guide provides only general suggestions regarding the geometry of the mesh to be used and no specific information on the density. According to the Moldflow Manual [7] , a 3D mesh requires a longer computational time, could provide better accuracy and is preferable for parts with thickness variation or mass accumulation (chunky regions). No further details regarding the most suitable mesh and density are given.
Miranda and Nogueira [8] studied the influence of gas entrapment on the plastic injection process using CAE and experimental analysis. Four thermocouples were used to measure the real mold wall temperatures with variations in the injection temperatures and pressure. For the best combination of parameters with an appropriate venting outlet the cycle time was 35% lower than without and there was a good fit between the experimental and simulation temperature data. The authors evaluated the mesh refinement and found that the smaller the element size the more reliable the geometry of the simulation domain will be.
The results of the simulation can be influenced by the quality of the mesh (in general, by the homogeneity of the distribution of the mesh elements) and its density (distance between two adjacent nodes). Both of these factors can affect the computation time and the precision of the simulation [9] . According to Kovács and Sikló [10] , a high number of nodes can cause a fluctuation in the results, probably due to truncation and rounding during the calculation. On the other hand, the use of a low-density mesh can lead to a deviation in the results since a low number of nodes does not represent adequately the geometry.
Although it cannot be considered a new technology, a scarcity of studies on the geometry and density of the mesh used to simulate the plastic injection molding process was verified during a review of the literature. Thus, the effects of these factors were investigated in this study. The aim of the research was to understand and quantify the influence of the mesh geometry and density in the CAE software used to simulate the plastic injection process. Critical variables inside the mold (pressure and temperature) were simulated and the results compared with experimental data.

Plastic Injection Molding Simulated Using CAE Software
The simulation process can be divided into three steps: i) pre-processing, where the mesh is generated and boundary conditions are applied; ii) processing, where the CAE software solve several equations in the mesh domain; and iii) post-processing: where the user interacts with the results, performing an evaluation and implementing changes if necessary. The main aspects of CAE software, the mesh generation (pre-processing) and the solving of equations (processing), will be discussed ahead.

Mesh generation
In the pre-processing stage, the CAE software (or other mesh generating software) generates a bidimensional or tridimensional mesh to describe the mold and/or the part domains of the plastic molding process [3] . In the case of a bidimensional mesh, the quality of the mesh is related to the distribution of the elements along the surface and the homogeneity on its length. The quality of a tetrahedral element is commonly measured in terms of minimal and maximal dihedral angles, since small or large angles result in less accuracy of the solution. The numerical solution stability is mainly affected by the worst tetrahedral element [11] .
Different approaches have been developed to generate a mesh as follows: a) The Delaunay approach attempts to distribute a set of vertices in the domain and additional vertices can be iteratively added as needed. In this method a 2D triangular mesh is usually applied, and it is not easily extended to the 3D complex domains [11,12] .
b) Octree-based methods subdivide the domain enclosing the given mesh recursively until a certain stopping criterion is reached [11] . In the modified-Octree technique, a CAD model is divided into various sizes of solid cubes and converted into tetrahedral elements [13] .
c) The advancing front technique (AFT) is based on the generation of a tetrahedron mesh from the front toward the interior, creating new elements until the entire solid geometry is filled with elements. In this approach, new nodes are inserted inside the domain to achieve the specified shapes and sizes for the mesh generation [11] . This is widely applied in 3D tetrahedral meshing and allows a better control of the elements generated [13] . Jin and Tanner [14] summarized the 3 main steps to run the AFT method, as follows: i) Discretize lines which form surfaces to generate the initial front.
ii) Discretize triangular faces on the surface to form the initial front.
iii) Discretize the triangular domain to create tetrahedral elements, only for elements on the advance front, checking the volume and upgrading mesh elements as necessary.
As described in Figure 1, to generate volumetric elements from a planar triangle it is necessary to create a new node (d) (offset from the face) using the existing nodes (a, b and c). Jin and Tanner [14] suggest limiting the new node location to within a sphere, with the center at d and a radius of 1.25δ, where δ is the edge reference dimension of the element to be created (Figure 1a). New nodes should not be too close to existing nodes, line segments or faces, and only new nodes whose distance from the other nodes is not less than 0.55δ can be selected. To avoid inconsistency in the calculation Jin and Tanner [14] sugestes that n abd . n abc > cos 75 and [15] that n abd . n abc > cos 60ᵒ, thus tetrahedron abcd is reasonably well-shaped ( Figure 1b).
Ito et al. [15] proposed a mathematical model to analyze the normal vector of the element and create the mesh as smooth as possible. To do so, the author used a Laplacian method with AFT to find the radius of all neighboring nodes to evaluate the triangles and correct them, in order to create a tetrahedral element which is as equilateral as possible. Parmar and Kaiser [16] evaluated the results for the temperature and the warpage obtained from simulations with an imported mesh compared with the use of 3D elements generated directly in commercial CAE software. The meshing method showed minimal variations and no influence on the simulation results. This study highlights the robustness of the mathematical mesh methods currently used in CAE software.

Solver equations
In this phase, using the mesh generated with proper boundary conditions, the software models the process and performs mathematical calculations. The main calculus on CAE software are the mold cooling (COOL), polymer injection (FILL) and packing (PACK), as well as dimensional calculations associated with shrinkage and warpage (WARP). In some commercial software programs, mechanical analysis of the mold and inserts can also be conducted.
For more than a decade, computer software has been available to solve the relevant equations for a Newtonian fluid flowing in a cold cavity. CAE software uses numerical methods, such as the finite element method (FEM) or finite difference method (FDM), to solve structural, thermal, fluid-dynamics and rheological equations. The predictions of the pressure curves inside the mold cavity can be significantly improved by introducing the effect of viscosity on pressure and the mold cavity deformation [17] . To solve thermal, stress and fluid-dynamics issues, the CAE software uses the elements domain, generated in the pre-processing step. To ensure accuracy, the boundary conditions applied must be close to experimental process conditions.
In the injection molding processing, the fluid is non-Newtonian with high viscosity (10 3 to 10 4 Pa.s) and the pressure is as high as 10 6 to 10 10 Pa [18] . In the generalized Newtonian model, the process can be described by mass, momentum and energy conservation equations. Although these equations describe different phenomenon, they come from the general transport equation (Equation 1): where ϕ and Γ are called characteristic variable and diffusivity, respectivally, and S ϕ is the source term [19] . The variables i u , j u and k u are the components of the velocity field in the three-dimensional Euclidian space with coordinates i x , j x and k x .
Zhou et al. [19] shown that by the appropriate choice of ϕ, Γ and S ϕ it is possible to obtain the conservative equations of the generalized Newtonian model. The alternatives to ϕ, Γ and S ϕ in the general transport equation are presented in the Table 1.
In Table 1, η and k are the viscosity and the thermal conductivity of the fluid respectively. The source terms for the momentum equations include the effects of pressure gradients and inertial forces. The source term of the energy equation is composed by the sum of the dissipated power due to viscous stress and other heat sources e S . The characteristic variable of the energy equation is the internal energy e of the fluid. However, it is common for this equation to be rewritten as a function of fluid temperature T. To deduce this expression, a relationship between internal energy, enthalpy h and measurable properties of the fluid is used. Zhou et al. [19] propose the thermodynamic relation: In energy equations the fluid is generally considered incompressible [19] . Considering this and the enthalpy as a linear function of temperature, we can obtain ( ) where p c is the specific heat at constant pressure. During COOL there is a zero velocity field, so the energy equation can be simplified to: The effective shear rate γ is a scalar variable derived from the shear rate tensor γ =  through the equation: where ( ) : is the double dot product.
The relation between the shear rate γ and the shear stress τ is given by a rheological equation. According to Fernandes et al. [20] , the viscosity ( ) is a function of the shear rate, the temperature and the pressure as show the Equation 6.
The rheological behavior has still not been fully defined and the equations are theoretically explained by the principles of continuum mechanics. However, extensive empirical data needs to be collected for the material characterization. In recent years, several models have been developed to describe the non-linear behavior of non-Newtonian liquids, such as polymers. The models used include the Cross-WLF Viscosity and Matrix viscosity model.
According to Li and Shen [21] , the calculations can be simplified by considering as incompressible Newtonian flow and neglecting the surface tension during FILL. During PACK the compressibility of the melt shall be taken into account. Therefore, a dependency model of the specific volume with temperature and pressure it's necessary. The modified Tait model is considered more suitable, because it can predict de abrupt volumetric change in the liquid-solid transition for amorphous polymers [20] . The modified Tait model is given by Equation 7.
where ( ) , v 0 T is the specific volume at zero gauge pressure,  [20] . According to Zhou et al. [19] simplifications can be made in the generalized Newtonian model to fit majority of injection cases, thin shell structures. Furthermore, due to the long molecular chain structure of polymers the viscous shear stress is much bigger than the inertial forces. As a result, velocity in the thickness direction and inertial forces can be neglected and the pressure is a function of planar coordinates [20] . Under these assumptions, we have the Hele-Shaw model: where: The boundary and initial conditions of the generalized Newtonian model and of the Hele-Shaw model include zero velocity on the normal and tangential directions of the mold cavity walls. Similarly, we have zero gradient pressure along the normal direction of the mold cavity walls. At the flow front and at the surface where melt enters the cavity (the gate) the pressure boundary conditions are p 0 = and inj p p = , respectivally. The injection pressure inj p and the gate speed are determined by the melt flow rate at the gate. The temperature boundary conditions can be prescribed in each boundary [19] .
To solve the boundary conditions problems described by the Hele-Shaw or the generalized Newtonian models, Moldflow's software uses the FEM. In terms of the general transport equation, the FEM approximates the unknown function ϕ over the domain of a finite element to obtain a system of algebraic equations. The approximation by a weighted procedure: and creates a residuo R: which should be zero at exact solution. The terms a N are known functions of local coordinates (basis functions) of the finite elements, while  a ϕ are unknown parameters. The best estimative of  a ϕ ensures that: where b W is an arbitrary function and n is the number of nodes in the finite element domain Ω . In order to avoid high-order derivatives, and achieve a weak formulation of the problem, integration by parts can be used [22] . The general transport equation is in the strong formulation due to the second-order derivative in the diffusion therm. W is the sum of a basis function b N and an artificial diffusion term that ensures the convergence for partial differential equations that not admit a weak formulation [22] .

Materials and Methods
To achieve the aims of this research, a mold was designed and manufactured. Temperature and pressure sensors were installed inside the mold cavity and a data acquisition system was developed to acquire the data (in real-time) during the injection molding process. Systematic CAE simulations were conducted in Moldflow Plastics Insight 2014 and analysis was carried out with the simulated results and the experimental data. The computation time was also recorded for all simulations.
The workpiece geometry designed for this study had a thickness of 2 mm, diameter of 140 mm and five equidistant cavities ( Figure 2). The cold sprue dimensions were: length of 60 mm, entrance with a diameter of 6.5 mm and draft angle of 2°. The cavity has a volume of 50.7 g/cm 3 . The cooling system had a diameter of 8 mm and was at a distance of 19 mm from the workpiece (U shape). The design of the mold was also aided by CAE.
Three CAE mesh geometries were evaluated: i) 2D midplane (MP); ii) 2D dual-domain (DD); and iii) 3D tetrahedral. The influence of the density of the mesh, which is the result of the maximum segment length defined by the user in the CAE software, was evaluated using four different values. The maximum segment length was evaluated with values from 1 mm to 6 mm using the same CAE tool tolerance of 0.1 mm for mesh generation. Table 2 shows the input variables and the number of elements generated in each mesh (density) according to the mesh geometry and maximum segment length.
On observing the extremes of the lengths evaluated for the 2D dual-domain and 3D meshes, it can be noted that the number of elements generated increases by around a factor of 7. Therefore, these two geometries are expected to have the greatest influence on the calculation time. The density for the 2D midplane geometry varied little as a function of the maximum segment length. Thus, using this mesh geometry, it is possible to obtain a higher precision without a large increase in the number of elements and, probably, little change in the computational time. Figure 3 shows an example of the mesh density for the 2D midplane geometry according to the maximum segment length.
The satisfactory quality of all meshes generated was verified by a computational tool available in the CAE software called 'Aspect Ratio', which calculates the level of distortion of the segment length of the meshes,   Figure 3. Overview of the mesh density according to the maximum segment length. and no significant distortions were found. The advancing front technique (AFT) was used to generate the 3D mesh. The simulation of the injection molding process was carried out in CAE software Moldflow Plastics Insight 2014, run on an up-to-date personal computer. In the simulation was selected the Cross-WLF model and the mold expansion was taken into account. A preliminary analysis to determine the COOL, FILL, PACK and WARP values was then conducted. The mold transient temperature regime was also evaluated for the 3D mesh. The experimental molding of the workpieces was conducted in a HAITIAN SA1200/410 machine. The switchover was setted for 98% of the volume filled in the simulations for the mold design, in the experiment short shoting technique was carried out and checked by the volume filled. The injection molding process parameters are shown in Table 3.
A piezoelectric sensor (Kistler, model 6190CA) was installed inside the mold cavity to acquire the pressure and temperature signals in real-time during the injection cycles ( Figure 4a). This sensor has a T-type thermocouple to measure the temperature of the melted material (in contact). The sensor was positioned 22 mm from the feed channel as shown in Figure 4b.
The sensor was calibrated for this application and the signals were amplified by a Kistler Amplifier 5039A221 and captured by an Agilent 34970 board.
The results for the CAE computation time and the temperature and pressure data were analyzed.

Computation time
The time required to compute the cases investigated can be seen in Figure 5 . As expected, the computation time increases with the mesh density (more elements). In the case of the 2D midplane mesh, using a maximum segment length of 1 mm, the simulation did not converge. A truncation might have led to an infinite calculation loop. Figure 5 also shows that the mesh density has a stronger influence on the computational time in the case of 3D elements compared with the 2D counterpart.
CAE users usually reduce the mesh density to shorten the computation time. However, in contrast to practical use, the results of Figure 5 show that, for all three mesh geometries, on using a mesh with over 2 mm of maximum segment length the computation time did not reduce significantly.
In general, the computation time for the 3D mesh was much longer (by up to a factor of 6) compared with the other meshes. However, the number of elements was around 10 times higher. This means that the calculation time for each element must be shorter for the 3D mesh. Endorsing this fact, varying the maximum segment length from 4 mm to 2 mm, the computation time increased by 21% and the number of elements increased by 43%. The accuracy of the simulated values according to the mesh are discussed below. Figure 6 shows the profiles for the dynamic pressure during the practical injection cycle and the simulated results. Firstly, it can be noted that for all cases the values for the simulated pressures increase before the experimental (real) values. This could be because of a delay in the drive of the injection machine (that could be up to 0.1 seconds according machine datasheet), when the control on the machine changes from the volumetric (speed) to pressure domain. This could not be properly considered in the simulation.

Analysis of the pressure data
The machine switchover delay, as can be observed in the Figure 6, resulted in a decrease of the preasure after the point 2 when compared to the simulation. This could also resulted in thiker frozen layers what explain the bigger preassure peaks comparing to midplane and dual-domain mesh simulation. Besides, the PP crystallization it's very sensible and can influence the pressure inside the cavity during cooling phase.
It is also observed that the simulated pressure drops faster than the experimental pressure, that is, in about half the time. This can be attributed to the viscoelastic behavior of the injected material. At the end of the packing phase, the material is subjected to an abrupt pressure variation followed by a relatively long relaxation period during the cooling phase, and the elastic component results in a delay to responses such as a drop in pressure. The simulation was not able to identify this phenomenon, which could affect the cycle time and the product quality.
In Figure 6, two important moments of the injection process can be identified: the maximum pressure point [7] and the peak at the end of the injection [13] . These specific moments are used to evaluate the pressure. Table 4 reports the simulated results for the pressure, the experimental pressure and the percentage of error between them.
The results show that there is no direct relation between the mesh geometry and density with the precision of simulation, which was an unexpected finding. Highlighting this situation, the most sophisticated mesh (3D) using the shortest maximum segment length (1 mm) generated the highest number of elements, but it resulted in the highest error (117%). Also, notably, this mesh required the longest computation time (around 6 h). The 2D dual-domain mesh results had the lowest deviations from the experimental data, with a minimum error of -3% and maximum of -12% (for the maximum pressure).
The 2D midplane mesh was associated with the shortest computation time and the lowest deviation of the experimental pressure (12% to 35%). This mesh calculated some peaks in the pressure at the end of the injection time (point 2) that were not present in the experimental injection molding. The computation for the 2D midplane with a maximum segment length of 1 mm did not converge. On analyzing the data error report of the simulation, this appears to be due to rounding and truncation, as suggested in the literature [10] .
The greater errors observed for the 2D midplane meshes compared to the 2D dual-domain meshes can be attributed to the behavior of the thin-walled injection flow. The melt can  be regarded as a general Hele-Shaw flow, which neglects the "fountain flow" phenomenon at the front of the melt [24] .

Analysis of the temperature data
Since the 2D meshes (midplane and dual-domain) are not able to simulate the dynamic temperature during the cycle, the temperature analysis was split into two parts:  i) static analysis at the end of the cycle; and ii) dynamic temperature during the cycle using a 3D mesh. In order to obtain a steady regime of the mold temperature, 10 batches of injection moldings were required. Figure 7 shows the results of the static analysis of the temperature.
The results obtained with the 2D dual-domain mesh did not vary with the mesh density and the error value was higher. The 2D midplane mesh was more accurate in this case, but the minimal error was around 14%, compared with the experimental (real) temperature. Figure 8 shows the dynamic temperature, obtained using the 3D mesh, during three injection cycles. It is possible to identify the start of the cycle, the end of the filling of the plastic material and the opening time.
In the case of the 3D mesh no significant influence of the maximum segment length was observed and the error between the simulated and the experimental temperature was lower than 2%. Therefore, the use of a 3D mesh with a longer segment length can save considerable computational time without loss of accuracy in the temperature simulation.
Another point observed in Figure 8 is that the simulated temperature rises sharply, in contrast to the experimental values. An angle α can be observed between the line of the simulated temperature and the line formed with the data obtained with the sensor, representing the divergence in the values. With the simulation is possible to evaluate, independently, the plastic mesh domain and mold mesh domain at the same instant of time. Thus, is possible to evaluate the effects of the thermal exchange on the temperature at specific points, where the mold and the plastic touch each other, as shown in Figure 9.
For the simulation results for the plastic and mold domain profiles, the temperature increase forms a vertical line, but differing in magnitude. The angle α represents the divergence between the simulation and experimental curves for the temperature. The resistance between the plastic material and the mold, and between the frozen layers of the polymer (results of the fountain flow phenomenon), does not allow instantaneous changes in temperature. The divergence of the simulated results from the sensor data can induce an error or a lack of convergence in a closed loop control, resulting in the need for manual adjustment of the parameters. The angle α identified in this study will be of use in further investigations. It should be noted that this divergence was observed only for the temperature simulation.

Conclusions
This paper reports the computational time and the temperature and pressure results obtained in simulations carried out with different CAE meshes. Experimental data was also acquired in real time. The main conclusions of this study can be detailed as follows: 1. The results show that the computational time was mainly influenced by the geometry of the mesh rather than the mesh density and it varied by up to a factor of 1,000 depending on the mesh geometry and density. In the case investigated, the mesh with the highest density (1 mm maximum segment length) was associated with the longest computational time.
2. In general, the use of a lower density mesh did not reduced significantly the accuracy of the simulation, which was unexpected, and the computation time was shorter.
3. The CAE simulation can reach a relatively good accuracy for the temperature parameter, when compared with the temperature measured experimentally inside the mold. For the 3D tetrahedral mesh the divergence was less than 2%. Therefore, in cases where high accuracy is required, the use of a 3D mesh is recommended.
4. Although the simulated temperature results showed good accuracy, the shape of the curves presented some divergence. The simulated temperature increased sharply, along a vertical line, during the injection period. However, in the real process the temperature increases in a nonlinear manner, as expected. This divergence was identified as angle α. This angle can affect the accuracy of the simulation and vary according to the material and the process parameters.
5. For the pressure parameter, the 3D mesh results show a significant divergence (of up to 117%) from the measured data. The closest results were obtained for the pressure simulated with the 2D dual-domain mesh, with a maximum error of 3%.
6. The pressure curves generated by the simulations were accurate up to the end of the speed control domain of the machine. The switchover to the pressure control domain resulted in errors in the simulations. The software was not able to properly identify the decompression during the process, indicating around half of the actual time required until the end of the packing phase, which can affect the filling and the integrity of the plastic part.
7. To have a more accurate simulation of the injection molding process by CAE systems, new developments should include the dynamic limitations of the process and the machine, such as the response time for temperature alterations, and input it into the software database.
To summarize the conclusions, for the pressure parameter, the 2D dual domain resulted in the most precise simulation and it was not influenced by the mesh density. Therefore, the use of this type of mesh and a lower density can provide a faster and more precise simulation of the pressure inside the mold. For the temperature parameter, the 3D mesh was more precise and it was also not influenced by the mesh density. Therefore, the use of a 3D mesh together with a lower density can allow a faster and more precise simulation of the temperature inside the mold.