Pressure Analysis of Dynamic Injection Molding and Process Parameter Optimization for Reducing Warpage of Injection Molded Products

Plastic injection molding technology is one of the important technologies for the manufacturing industry. In this paper, a numerical dynamic injection molding technology (DIMT) is presented, which is based on the finite element (FE) method. This numerical simulation method introduces a vibrational force into the conventional injection molding technology (CIMT). Some meaningful work has been executed for investigating the mechanical evolution behavior of DIMT. As the basis for illustrating the mechanism in warpage optimization results, dynamic parameter analysis on the rule of pressure response is performed in detail. In the warpage optimization work, three kinds of structure with different molding materials are selected as the comparison. The final warpage of each product is efficiently minimized by using a Gaussian process-based sequential optimization method. From the further discussions, the features of DIMT in improving the molding quality are revealed, indicating that it may not be appropriate for geometrically large structures. This study has significant meaning for the actual injection molding industry.


Introduction
Plastic injection molding technology has been regarded as one of the most important technologies for rapid, efficient manufacturing. During the injection molding process, a plastic melt is injected into the void cavity by the screw of the injection molding machine, and subsequently solidified in the cavity after a certain cooling time. It is a complicated multi-physical process involving a thermal and mechanical evolution history. The level of molding defects, like warpage or shrinkage magnitude, length of weld line, etc., are associated with this complicated thermal and mechanical history. Therefore, it has attracted researchers' attention to investigate the thermal and mechanical behavior in the injection molding process, and relevant process optimization methods to improve the molding quality for a long time.
In fact, the thermal and mechanical behavior of the plastic material depends on the process parameters, molding temperature, pressure, etc. As a special injection molding technology, dynamic where ρ, c p , β, and k represent four physical properties, density, specific heat, coefficient of thermal expansion, and thermal conductivity, respectively. v and f represent the velocity field and body force per unit mass, respectively. σ is the stress tensor that can be decomposed into an isotropic part and a non-isotropic part, as in Equation (4).
. q is the heat generation rate per unit mass. κ is the expansion viscosity that can determine the heat generation level during expansion. T and p e represent the temperature and equilibrium pressure, respectively.
Obviously, it is difficult to solve the complete equations above, so we have to make some physical assumptions for simulating the plastic injection molding process. The main assumptions are: (a) For thin-walled injection molded products, the flow behavior of the melt can be regarded as the general Hele-Shaw flow, which neglects the "fountain flow" phenomenon at the front of the melt. The finite mesh type is the mid-plane mesh with triangular elements. (b) Considering the high viscosity effect of the plastic melt, inertia force and gravity are neglected. Therefore, the flow process of the melt is treated as a laminar flow, and regarded as the general Newtonian fluid. The melt is dominated by shear stress. (c) No-slip assumption. This means that there is no relevant movement between the melt and the mold cavity surface, or between the melt and the solid plastic material surface. (d) For the energy transport issue, there is no heat source in the melt in practice. We only consider heat conduction, but neglect heat convection along the thickness direction. In contrast, we consider heat convection and neglect heat conduction along the flow direction.
According to Assumption (b), the constitutive equation of the melt should coincide with the form of the Newtonian fluid, which reads as: where . γ ij is the rate of strain tensor. η(T, p, . γ) is the viscosity of the melt, and is a highly non-linear function on temperature T, pressure p, and shear rate . γ, which is defined as: Polymers 2017, 9, 85 4 of 23 Many researchers made efforts in order to obtain a perfect model of η(T, p, . γ) [21][22][23], listed in Table 1. However, the Cross-WLF model is suitable to be used as a computation model for most plastic materials. In this paper, Cross-WLF model is employed. Table 1. Viscosity models for describing the non-Newtonian behavior of plastic melt [21][22][23].

Model
Contributor By means of integrating the continuous Equation (1) along with the thickness and Galerkin-weighted residual method, we can obtain Equation (13) completely for solving the flow process: where: where ξ is the position between the melt and solid material along the thickness.
In addition, many properties of the plastic material are variable with temperature and pressure, like thermal conductivity and density. For describing the thermal conductivity evolution, we employ: where λ 1 to λ 5 are constant parameters. The pressure-volume-temperature relationship, which is usually called the pvT state equation, is the critical rule to characterize the change of density (specific volume) with respect to temperature and pressure. For the pvT state equation, there are several state models to describe the property change, like the Spencer-Gilmore (S-G) model [24] and the two-domain Tait model [25]. In this study, we employ the two-domain Tait model, as Equation (15) expresses: where C = 0.0894, and: where b 1 to b 9 are material parameters, T t represents the transition temperature, which is a crystalline temperature for crystalline plastic and the glass transition temperature for amorphous plastic. As an explanation of the pvT relationship, Figure 1 presents the specific volume change of PC Lexan 105 with respect to temperature and pressure.
Polymers 2017, 9, 85 6 of 25 explanation of the pvT relationship, Figure 1 presents the specific volume change of PC Lexan 105 with respect to temperature and pressure.

The Finite Element Equation for Computing the Pressure Field
For simulating the flow process during injecting, we need to deduce the finite element discrete form of Equation (13). Consider that the pressure in any element is interpolated by: where index m indicates the node number of the current element. Lm is the interpolation function,

The Finite Element Equation for Computing the Pressure Field
For simulating the flow process during injecting, we need to deduce the finite element discrete form of Equation (13). Consider that the pressure in any element is interpolated by: (20) where index m indicates the node number of the current element. L m is the interpolation function, which satisfies: where A is the area of the triangle element. Then, we can obtain the finial Equation (22): where K, P, and Q represent the stiffness matrix, pressure vector, and the net flow rate vector, respectively. Subscript G, M, and F represent the injecting gate, melt region and the flow front of the melt, respectively. Therefore, Q M = 0 in the melt, and P F = 0 at the front of the melt. For the conventional injection molding, flow rate Q G is usually a given constant during filling analysis, and P G is a given constant or linear form during the packing phase. We can make a slight change of Equation (22): for computing the pressure and velocity field in filling analysis with Q G given, and: for computing the pressure and velocity field in packing analysis with P G given.
For the dynamic filling phase, the flow rate can be deduced by the advancing displacement of the screw which is: where v ave is the average advancing velocity of the screw, S amp , f fill , and φ fill are amplitude, frequency, and phase of a vibration filling. Then, the flow rate is: where V is the volume of the product, A sc here is the area of the screw, t fill is the approximate filling time. ς equals (S amp ·2π·f fill ·t fill ·A sc )/V herein.
Although Equation (26) is available to be the dynamic filling process, there is a limitation that, in order to guarantee no melt reflux problems, (S amp ·2π·f fill ·t fill ·A sc )/V ≤ 1 has to be satisfied. This indicates that high frequency or high amplitude vibration cannot be realized, since if S amp or f fill is higher, f fill or S amp will be so low that the vibration will be "no sense" in practice. Therefore, we prefer to choose another case of dynamic filling rather than Equation (26) in the following study.
Therefore, if we express the advancing velocity of the screw in a simple form as: where ε is called the amplitude factor, we can easily obtain the dynamic filling flow rate expression as: Equation (28) avoids the constraint problem mentioned above, and can be realized easily as well. For the packing phase, making the similar manipulation the filling phase, the packing pressure curve is controlled by: where P final represents the finial pressure during the filling phase, P start is the packing pressure fraction at the beginning, δ is a factor to control the trend of the packing curve, θ, f pack , and φ pack represent the amplitude, frequency, and phase of a vibration term, respectively. Overall, with Equations (28) and (29), Equations (23) and (24) can be solved to investigate the pressure response analysis and the molding quality improvement analysis of plastic product of numerical DIMT, which will be argued in the following sections.

Pressure Response Analysis in the Cavity
Defects, like warpage, shrinkage, etc., depend on the mechanical and thermal evolution history. By integrating Equation (7) associated with the solution conditions, we can deduce the velocity formula as: where h is the thickness of the product. Obviously, solving the velocity field requires us to firstly solve the pressure field. Moreover, computing the heat contribution from shear motion relies on the pressure field as well. Therefore, pressure in the melt is a very important mechanical quantity to lead the material motion and the thermal evolution. This indicates that pressure response analysis can help us to understand the mechanical and thermal evolving history.
Therefore, a finite element model of a low-density polyethylene (LDPE) circle disk product (depicted as Figure 2) is modeled as the investigating objective, which is the same as the product of Yin's experiment [7]. The basic material properties of LDPE are listed as Table 2. The injecting gate is located at the center of the disk, which is colored with pink in the finite element model. mechanical and thermal evolving history.
Therefore, a finite element model of a low-density polyethylene (LDPE) circle disk product (depicted as Figure 2) is modeled as the investigating objective, which is the same as the product of Yin's experiment [7]. The basic material properties of LDPE are listed as Table 2. The injecting gate is located at the center of the disk, which is colored with pink in the finite element model.  In order to compare the simulation result of DIMT with the experimental observation [7], three monitor points, except the gate, A, B, and C, are selected in Figure 2b. The distance of OA, OB, and OC are 10.08, 29.40, and 46.34 mm, respectively. The melt temperature and mold temperature are 190 °C and 50 °C , respectively. Moreover, injecting time and packing time are both 2 s. Other dynamic parameters are listed in Table 3.  In order to compare the simulation result of DIMT with the experimental observation [7], three monitor points, except the gate, A, B, and C, are selected in Figure 2b. The distance of OA, OB, and OC are 10.08, 29.40, and 46.34 mm, respectively. The melt temperature and mold temperature are 190 • C and 50 • C, respectively. Moreover, injecting time and packing time are both 2 s. Other dynamic parameters are listed in Table 3. By finite computation, we can obtain the pressure response of these four monitors both for the CIMT case and DIMT case. From Figures 3 and 4, it is clarified that, compared with the experimental observation, our pressure results are efficient, although there is a certain difference between the computation and experiment. The value of the simulation results is slightly lower than the experimental value. This difference is mainly because the material property and processing parameters of computation cannot be exactly the same as the actual experiment. However, this indicates that our computation method can be used as a simulation tool for investigating the mechanical evolution and the product quality under DIMT further. observation, our pressure results are efficient, although there is a certain difference between the computation and experiment. The value of the simulation results is slightly lower than the experimental value. This difference is mainly because the material property and processing parameters of computation cannot be exactly the same as the actual experiment. However, this indicates that our computation method can be used as a simulation tool for investigating the mechanical evolution and the product quality under DIMT further.   observation, our pressure results are efficient, although there is a certain difference between the computation and experiment. The value of the simulation results is slightly lower than the experimental value. This difference is mainly because the material property and processing parameters of computation cannot be exactly the same as the actual experiment. However, this indicates that our computation method can be used as a simulation tool for investigating the mechanical evolution and the product quality under DIMT further.   There are several interesting phenomena in Figure 4. The first one is that the pressure response in the cavity has a fluctuation compared with the CIMT case. The second one is that this fluctuation amplitude gets larger and larger as time goes by for each monitor. The third one is that, although the observation of the pressure response rule for different monitors is the same, the further the monitor location is away from the gate, the smaller the fluctuation amplitude is. The fourth one is that, if we observe the fluctuation frequency, we find that this frequency seems like the one of the dynamic flow rate. If we make a fast Fourier transform (FFT) operation to this fluctuation, we can obtain the frequency domain as Figure 5, and the result shows that this fluctuation frequency is definitely the same as that of the flow rate. These observations exist in the filling phase. In the packing phase, the mean difference between these monitors is the fluctuation amplitude.
Recalling the finite element equation for computing the mechanical field under DIMT Equation (23), we can find the explanation of all of these phenomena. Actually, the entire solution of the equations not only includes the solution of the CIMT part, but also superposes the solution of the "vibration term". For the solution of the CIMT part, it gives people the main trend of the pressure response curve, whereas the solution of the "vibration term" makes a certain pressure revision on the basis of the CIMT solution. For example, due to the periodic behavior of the flow rate in the dynamic filling phase, the higher flow rate of the DIMT case than the constant one of the CIMT case induces a pressure ascent, while the lower flow rate induces a pressure descent. This is the main reason that the pressure of the DIMT case fluctuates around the pressure of the CIMT case. For the difference of the fluctuation amplitude between these monitors, it is because the pressure value depends on the flow length between the node and the flow front. The larger the flow length is, the larger the fluctuation amplitude presents.
equations not only includes the solution of the CIMT part, but also superposes the solution of the "vibration term". For the solution of the CIMT part, it gives people the main trend of the pressure response curve, whereas the solution of the "vibration term" makes a certain pressure revision on the basis of the CIMT solution. For example, due to the periodic behavior of the flow rate in the dynamic filling phase, the higher flow rate of the DIMT case than the constant one of the CIMT case induces a pressure ascent, while the lower flow rate induces a pressure descent. This is the main reason that the pressure of the DIMT case fluctuates around the pressure of the CIMT case. For the difference of the fluctuation amplitude between these monitors, it is because the pressure value depends on the flow length between the node and the flow front. The larger the flow length is, the larger the fluctuation amplitude presents.

Pressure Response Analysis at the Gate
Pressure fluctuation exists in the cavity, while it is the most obvious at the gate. Therefore, we can focus on the injecting pressure to investigate the influence of DIMT parameters on the evolution of the pressure field.

Influence of the Filling Amplitude on the Transient Injecting Pressure
In this subsection, we make ffill equal to 6 Hz, φfill equals 0, then analyze that the transient injecting pressure response of nine different values for ε from 0 to 0.8 at 0.1 intervals. The injecting pressure curves are depicted as Figure 6. This shows that, as ε increases from the conventional case, the volatility level increases. This observation also exists for other cases of ffill as shown in Figure 7. This is because the pressure solution of DIMT is the revision of CIMT case from Equation (23), which

Pressure Response Analysis at the Gate
Pressure fluctuation exists in the cavity, while it is the most obvious at the gate. Therefore, we can focus on the injecting pressure to investigate the influence of DIMT parameters on the evolution of the pressure field.

Influence of the Filling Amplitude on the Transient Injecting Pressure
In this subsection, we make f fill equal to 6 Hz, φ fill equals 0, then analyze that the transient injecting pressure response of nine different values for ε from 0 to 0.8 at 0.1 intervals. The injecting pressure curves are depicted as Figure 6. This shows that, as ε increases from the conventional case, the volatility level increases. This observation also exists for other cases of f fill as shown in Figure 7. This is because the pressure solution of DIMT is the revision of CIMT case from Equation (23), which was explained in Section 3. However, it also implies from the figure that if ε gets closer to 0, the "dynamic effect" will become inefficient, since the dynamic flow rate yields into the conventional case, and parameters f fill and φ fill become useless process parameters simultaneously.
Polymers 2017, 9,85 11 of 25 was explained in Section 3. However, it also implies from the figure that if ε gets closer to 0, the "dynamic effect" will become inefficient, since the dynamic flow rate yields into the conventional case, and parameters ffill and φfill become useless process parameters simultaneously.

Influence of the Filling Frequency on the Transient Injecting Pressure
In this subsection, we choose ε = 0.5, φfill = 0, ffill belongs to [0, 12] Hz of 2 Hz interval, then, we simulate the filling process for each case of seven different values of ffill. The final transient pressure curves of the front four cases are collected in Figure 8. According to the analysis in Section 3, we know that ffill actually affects the frequency of the dynamic injecting pressure by using FFT technology. In Figure 8, it is notable that for the given ε, the entire volatility level of injecting

Influence of the Filling Frequency on the Transient Injecting Pressure
In this subsection, we choose ε = 0.5, φ fill = 0, f fill belongs to [0, 12] Hz of 2 Hz interval, then, we simulate the filling process for each case of seven different values of f fill . The final transient pressure curves of the front four cases are collected in Figure 8. According to the analysis in Section 3, we know that f fill actually affects the frequency of the dynamic injecting pressure by using FFT technology. In Figure 8, it is notable that for the given ε, the entire volatility level of injecting pressure for each case is almost the same, even though f fill in each case is different. However, different f fill can cause different injecting pressures at the end of the filling phase. This is not a good expectation in some cases. For instance, if ε is close to 1 and f fill takes a small value, it will cause the short shot issue due to low injecting pressure when the flow rate descends below the mean value of the CIMT case at the late period of filling. pressure for each case is almost the same, even though ffill in each case is different. However, different ffill can cause different injecting pressures at the end of the filling phase. This is not a good expectation in some cases. For instance, if ε is close to 1 and ffill takes a small value, it will cause the short shot issue due to low injecting pressure when the flow rate descends below the mean value of the CIMT case at the late period of filling.

Influence of the DIMT Parameters on the Mean Injecting Pressure
From the previous observations, we know that pressure of DIMT fluctuates around the pressure curve of CIMT, sometimes higher and sometimes lower. It is not valuable to compare the pressure between DIMT and CIMT at any given time, but valuable to study the influence of dynamic parameters on the mean pressure during filling phase. The mean pressure during the filling phase

Influence of the DIMT Parameters on the Mean Injecting Pressure
From the previous observations, we know that pressure of DIMT fluctuates around the pressure curve of CIMT, sometimes higher and sometimes lower. It is not valuable to compare the pressure between DIMT and CIMT at any given time, but valuable to study the influence of dynamic parameters on the mean pressure during filling phase. The mean pressure during the filling phase represents the fluidity of the melt and the total input work, which can be defined as: By using Equation (32), we have computed the mean pressure of the different filling frequencies, and obtained the relationship between the mean pressure and ε, which is shown in Figure 9. Obviously, the mean pressure of the DIMT case is lower than that of the CIMT case from the figure. Moreover, as the amplitude increases, the mean pressure decreases for each one of the investigating frequency simulations. This indicates that the melt flows much easier under DIMT than CIMT, and more energy is saved for manufacturing the product. This is good news for producers since the energy saved is much more considerable for continuously manufacturing products. However, this effect is not obvious for ε belongs to [0, 0.3] since small ε induces an unimpressive fluctuation, as mentioned in Section 3. If we change the upper limit of integral tfill into t, which belongs to [0, tfill] in Equation (32), we can observe the evolution of the mean pressure. For instance, Figure 10 presents the mean pressure curves for ffill = 6 Hz. In this figure, it reveals that the conclusion in Figure 9 of the above study cannot always hold during the entire filling process, but holds after a certain time in the filling process. At the beginning of the filling phase, there is a wing-shaped region which is enclosed by the mean pressure curve of DIMT and that of CIMT. The area of this region increases as ε increases. Moreover, for a given ε, the mean pressure curve presents a fluctuation. These phenomena also exist in other frequency cases, as shown in Figure 11. If we change the upper limit of integral t fill into t, which belongs to [0, t fill ] in Equation (32), we can observe the evolution of the mean pressure. For instance, Figure 10 presents the mean pressure curves for f fill = 6 Hz. In this figure, it reveals that the conclusion in Figure 9 of the above study cannot always hold during the entire filling process, but holds after a certain time in the filling process. At the beginning of the filling phase, there is a wing-shaped region which is enclosed by the mean pressure curve of DIMT and that of CIMT. The area of this region increases as ε increases. Moreover, for a given ε, the mean pressure curve presents a fluctuation. These phenomena also exist in other frequency cases, as shown in Figure 11. pressure curve of DIMT and that of CIMT. The area of this region increases as ε increases. Moreover, for a given ε, the mean pressure curve presents a fluctuation. These phenomena also exist in other frequency cases, as shown in Figure 11.  At the beginning of the filling phase, it shows that the flow rate is increasing and is larger than the conventional flow rate from Equation (28) for this study. Therefore, the greater ε can cause a greater difference between the pressure curves of DIMT and CIMT, which makes a wing-shaped region. As time goes on, the effect of DIMT emerges that reduces the flow resistance, and the mean pressure descends below that of CIMT.
In order to explain the mean pressure fluctuation observation, we need to analyze the mathematical monotony of the mean pressure curve, which is deduced from Equation (33) as: At the beginning of the filling phase, it shows that the flow rate is increasing and is larger than the conventional flow rate from Equation (28) for this study. Therefore, the greater ε can cause a greater difference between the pressure curves of DIMT and CIMT, which makes a wing-shaped region. As time goes on, the effect of DIMT emerges that reduces the flow resistance, and the mean pressure descends below that of CIMT.
In order to explain the mean pressure fluctuation observation, we need to analyze the mathematical monotony of the mean pressure curve, which is deduced from Equation (33) as: where F(t) = t · p(t) − t 0 p(t )dt determines the monotony of the mean pressure curve. By using the mean value theory of integral, which states: where τ belongs to [0, t), and τ = 0 only if t = 0 for this problem. Then, Equation (34) is simplified as: Equation (35) associates the monotony of the mean pressure with the injecting pressure value. For CIMT, it is constant practice that p(τ) < p(t) since the injecting pressure is monotone increasing for filling the cavity. However, for DIMT, it does not hold anymore. Recalling the observation of the injecting pressure of DIMT, p(τ) may be smaller than p(t), and p(τ) may be greater than p(t), as well, due to the pressure fluctuation of DIMT. Therefore, the monotony of the mean pressure for DIMT is variable, which induces the fluctuation phenomenon for all cases in Figure 11.
It can be clarified from another point of view that if computing the gradient of F(t) with respect to t, we can obtain: There is no doubt that F(0) = 0. Then, Equation (34) can be rewritten as: Using the mean value theory of integral again, we can deduce a much simpler formula: Equation (39) associates the monotony of the mean pressure with the monotony of the injecting pressure. It can illustrate the phenomena in Figure 11 more directly with the help of the observations on the injecting pressure response.

Optimization Model
For thin-walled products, warpage is one of the most important molding defects. In the following sections, the purpose is to minimize the warpage of thin-walled products and argue the actual effect of DIMT compared with CIMT. The warpage optimization model is proposed as: where d i x (x), d i y (x), and d i z (x) represent the deformation in the x, y, and z direction of node i. N node and N var represent the node number and the variable number, respectively.
In order to compare the improvement level of warpage, the process parameter optimizations for CIMT and DIMT are both executed. For CIMT, melt temperature x 1 (T melt ), mold temperature x 2 (T mold ), filling time x 3 (t fill ), packing time x 4 (t pack ), initial packing fraction x 5 (P start ), and trend coefficient x 6 (δ) are selected as the optimization parameters. For DIMT, except the previous six parameters, filling amplitude x 7 (ε), filling frequency x 8 (f fill ), filling phase x 9 (φ fill ), packing amplitude x 10 (θ), packing frequency x 11 (f pack ), and packing phase x 12 (φ pack ) are chosen, additionally.
The optimization objective in Equation (40) is a "black" function with respect to x. In this paper, a Kriging surrogate model [26,27] is taken into consideration to evaluate the objective approximately, since it can give us more prediction information, and a maximum experiment improvement (MEI) function-based sequential optimization method [28] is employed to solve Equation (40). The detailed contents of the optimization method and procedure can be found in Appendix A. Then, the equal optimization issue changes to be: with the convergence criterion: where y I max and y I min are the maximum and minimum response values in the sample set of the I th searching for Equation (41), respectively. Herein, o c is 0.001 in this paper.

Structure and Parameters
In this section, three types of warpage optimization issues are studied to argue the improvement capability of DIMT compared with CIMT.

(a) Small Size Product with No Holes
The first product model is a polycarbonate (PC, the type is X-1) light cover with no holes in it. The geometry and the finite element model of the product are sketched as shown in Figure 12. The basic properties of the molding material are listed in in detail Table 4. The value ranges of 12 process parameters are tabulated as Table 5. For the optimization of CIMT, the design variables involved are the first six variables in Table 5. In order to compare the improvement level of warpage, the process parameter optimizations for CIMT and DIMT are both executed. For CIMT, melt temperature x 1 (Tmelt), mold temperature x 2 (Tmold), filling time x 3 (tfill), packing time x 4 (tpack), initial packing fraction x 5 (Pstart), and trend coefficient x 6 (δ) are selected as the optimization parameters. For DIMT, except the previous six parameters, filling amplitude x 7 (ε), filling frequency x 8 (ffill), filling phase x 9 (φfill), packing amplitude x 10 (θ), packing frequency x 11 (fpack), and packing phase x 12 (φpack) are chosen, additionally.
The optimization objective in Equation (40) is a "black" function with respect to x. In this paper, a Kriging surrogate model [26,27] is taken into consideration to evaluate the objective approximately, since it can give us more prediction information, and a maximum experiment improvement (MEI) function-based sequential optimization method [28] is employed to solve Equation (40). The detailed contents of the optimization method and procedure can be found in Appendix A. Then, the equal optimization issue changes to be:

Structure and Parameters
In this section, three types of warpage optimization issues are studied to argue the improvement capability of DIMT compared with CIMT.

(a) Small Size Product with No Holes
The first product model is a polycarbonate (PC, the type is X-1) light cover with no holes in it. The geometry and the finite element model of the product are sketched as shown in Figure 12. The basic properties of the molding material are listed in in detail Table 4. The value ranges of 12 process parameters are tabulated as Table 5. For the optimization of CIMT, the design variables involved are the first six variables in Table 5.    The second product model is a phone cover used in our previous work [20], which is made by the blend of PC and acrylonitrile butadiene styrene (ABS) copolymers. Herein, it is only combined as a comparison with other issues. In this case, there are multiple holes in the product compared with issue (a). Figure 13 is the geometry and the finite element model of the product. The basic properties of the molding material are listed in Table 6. The value ranges of 12 process parameters are tabulated as Table 7, in which the design variables involved in CIMT optimization issue are the first six variables.   The second product model is a phone cover used in our previous work [20], which is made by the blend of PC and acrylonitrile butadiene styrene (ABS) copolymers. Herein, it is only combined as a comparison with other issues. In this case, there are multiple holes in the product compared with issue (a). Figure 13 is the geometry and the finite element model of the product. The basic properties of the molding material are listed in Table 6. The value ranges of 12 process parameters are tabulated as Table 7, in which the design variables involved in CIMT optimization issue are the first six variables.   The third product is a large, narrow PC light cover. The geometrical structure information and finite element model of this PC cover is sketched in Figure 14. The basic properties information of   The third product is a large, narrow PC light cover. The geometrical structure information and finite element model of this PC cover is sketched in Figure 14. The basic properties information of the molding material is listed in Table 8. The value ranges of 12 process parameters are tabulated as Table 9. For the CIMT case, the design variables take the first six variables in Table 9.
Polymers 2017, 9,85 18 of 25 the molding material is listed in Table 8. The value ranges of 12 process parameters are tabulated as Table 9. For the CIMT case, the design variables take the first six variables in Table 9.

Results and Discussions
By using MEI-based optimization method in Appendix A, all of the issues of warpage optimization are solved. The optimization results of three CIMT cases are listed in Table 10, Table 12,  and Table 14, respectively. The optimization results of three DIMT cases are listed in Table 11, Table  13, and Table 15, respectively. For issue (a), the warpage of the PC light cover under DIMT (0.168 mm) is smaller than that of the CIMT case (0.197 mm), which is reduced by 14.72%. For issue (b), the warpage of the PC/ABS phone cover under DIMT (0.094 mm) is smaller than that of the CIMT case (0.221 mm), which is reduced by 57.47%. However, for issue (c), the warpage of the PC car light cover under DIMT (0.375 mm) is almost the same as that of the CIMT case (0.388 mm). This shows that DIMT do not improve the molding quality better than CIMT for this large narrow product. The warpage contour figures of all of the results are organized in Figures 15-17, respectively.

Results and Discussions
By using MEI-based optimization method in Appendix A, all of the issues of warpage optimization are solved. The optimization results of three CIMT cases are listed in Tables Tables 10, 12 and 14, respectively. The optimization results of three DIMT cases are listed in Tables Tables 11, 13 and 15, respectively. For issue (a), the warpage of the PC light cover under DIMT (0.168 mm) is smaller than that of the CIMT case (0.197 mm), which is reduced by 14.72%. For issue (b), the warpage of the PC/ABS phone cover under DIMT (0.094 mm) is smaller than that of the CIMT case (0.221 mm), which is reduced by 57.47%. However, for issue (c), the warpage of the PC car light cover under DIMT (0.375 mm) is almost the same as that of the CIMT case (0.388 mm). This shows that DIMT do not improve the molding quality better than CIMT for this large narrow product. The warpage contour figures of all of the results are organized in Figures 15-17, respectively.            For these three issues, the results show that DIMT can greatly reduce the warpage of the products, even though the structure and molding material are different. For issues (a) and (c), the dynamic parameters are all "efficient". However, the f pack and φ pack are both "inefficient" in issue (b) since the packing amplitude is 0, and the molding process yields into the CIMT case during packing phase. This is to say that, in issue (b), only dynamic filling parameters worked. The discussion of this optimization result is not the main work here, and has been stated in our previous work [20].
From the warpage contours, it shows that, in issues (a) and (b), the distribution area of the smaller warpage region under DIMT are larger than that of the CIMT cases, and the distribution area of the larger warpage region under DIMT is smaller than that of the CIMT cases. Therefore, DIMT can not only reduce warpage much more than CIMT, but also improve the entire molding quality for the small size products. Unfortunately, it does not happen for issue (c), which makes the molding quality much worse than CIMT. We should recall the pressure response observation of the cavity in Section 3. From the analysis in Section 3, we should realize that the pressure volatility relies on the flow length. If the geometry size is large, at the location far away from the gate, the pressure will not be affected by the vibration. This causes a great difference between the region around the gate and the region away from the gate. It may not have a good mechanical history for improving the molding quality. Therefore, it fails.
Actually, the in-cavity residual stress is the final result of the thermal and mechanical evolution before ejection. The uneven distribution of the in-cavity residual stress is the only reason for the warpage after ejection. DIMT takes an external force into the flow history. For small size product, it can disturb the flow behavior and may eliminate the magnitude and distribution difference of the in-cavity residual stress, which results in a lower warpage level compared with CIMT. Whereas, for large-sized, narrow products, this contribution is limited around the gate, and may not be good for improving the molding quality, even making it worse. Therefore, using DIMT should be seriously considered to produce the polymer product.

Conclusions
In this paper, the numerical DIMT is proposed by using the finite element method. More detailed work has been performed for observing and analyzing the mechanical evolution behavior of DIMT. The effect of dynamic parameters on the transient pressure of the cavity and the mean pressure of the gate are analyzed for realizing the mechanical behavior of DIMT. It is also fundamental for explaining the mechanism in the following warpage optimization work.
Three kinds of optimization issues with different molding materials are constructed in this paper. By using the EGO sequential optimization method, warpage values of all of the products are efficiently reduced. The further analysis and discussions reveal that DIMT is not available for all types of structures, especially for the large-scale, narrow structures. The main reason for this is that the thermal and mechanical evolution history can determine the in-cavity residual stress in the product before being ejected out of the mold. For small-sized products, dynamic process parameters make the key contribution to disturb the in-cavity residual stress distribution, which can finally increase the region of lower warpage and reduce the warpage difference in the product. Whereas, for large and narrow products, the dynamic contribution concentrates around the gate. It makes a great mechanical difference for this kind of structure, and causes a significant difference in warpage in the product. Therefore, as a recommendation, for this type of product, DIMT should be seriously considered. Overall, the content of this study focuses on revealing the mechanical mechanism of DIMT, and may be a guide for the actual manufacturing with DIMT.
Author Contributions: Xinyu Wang, Junfeng Gu, Changyu Shen proposed the idea of numerical DIMT, and developed the numerical method of DIMT, Xinyu Wang coded the procedure and wrote the paper, Hongxia Li, Shilun Ruan and Minjie Wang were responsible for analyzing the results. Zheng Li was responsible for the EGO method.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Kriging Surrogate and EI-Based Sequential Optimization Method (Equation Section (Next))
An optimization issue of the black function y(x) is given by: Find : x = x 1 , x 2 , x 3 , . . . , x m Min : y(x) S.T. : x i lb ≤ x i ≤ x i ub , x i ∈ x, i = 1, 2, . . . , m g j (x) ≤ 0, j = 1, 2, . . . , N g h k (x) = 0, k = 1, 2, . . . , N h where m is the number of the variable, and N g and N h are the numbers of the constraints.
It needs two datasets for constructing the kriging model for Equation (A1), the sample set X of the design variables and the corresponding objective value set Y. For the given variable sample set X and the corresponding objective set Y, there are two parts, a polynomial part and a stochastic function, in the kriging model shown in Equation (A2): where β k and f k (x) are the kth regression coefficient and polynomial. z(x) is the Gaussian stochastic function which follows norm [0, σ 2 ]. For x i and x j , z(x i ) and z(x j ) are correlated with the spatial weighted distance between x i and x j . The covariance is expressed as: Cov z(x i ), z(x j ) = σ 2 · R ij (λ, where R ij is a correlation function. We usually employ Gaussian form as shown in Equation (A4): where λ n is the nth correlated coefficient for the nth variables. The ordinary kriging (OK) model is commonly used in most studies, which f k (x) in Equation (A2) only employs the zero-order. The prediction and variance of the OK model at any point x* yields: where u = (l T R −1 l) −1 (l T R −1 Y). Y is the objective vector. r i = R(λ, x*, x i ). R is correlation function matrix. The element of l equals one. σ 2 is the maximum likelihood of σ 2 .
According to the properties of kriging model, we can define some optimal searching functions to exploiting the solution of problem Equation (40). The most efficient searching function is the expected improvement (EI) function, which is the expectation of the I(x) = y min − y(x) in the positive region of I(x). The detail formula of EI function is defined as: where y min is the minimum value in the current sample set Y.
In mathematics, E[I(x)] is an expectation surface to exploit a smaller design than the minimum in the current sample response set, or the most uncertain space design point from the expression Equation (A7). Therefore, Equation (A7) has both the local and global searching ability, and the maximum of E[I(x)] (MEI) contains the strongest potential to find the optimal solution, or points out the lack of information of the current surrogate model. This means that there exists little space for improvement to find a better design point than the current one, or to support a valuable spatial point to revise the current surrogate model. The MEI-based optimization method is usually called the efficient global optimization (EGO). Herein, we employ the MEI-based optimization method to solve the warpage optimization problem. Then the optimization problem changes into: Find : x = x 1 , x 2 , x 3 , . . . , x m Max : EI(x) S.T. : x i lb ≤ x i ≤ x i ub , x i ∈ x, i = 1, 2, . . . , m g j (x) ≤ 0, j = 1, 2, . . . , N g h k (x) = 0, k = 1, 2, . . . , N h (A8) The optimizing procedure is stated as: (a) Select the proper initial sample design X by using the mathematical sampling method, like orthogonal Latin hypercube sampling, optimized Latin hypercube sampling. etc., and compute the corresponding objective set Y; (b) Construct the kriging surrogate model with X and Y; (c) Construct the EI function and the optimization problem Equation (A8), and deal with the constraint conditions; (d) Find the optimal design x* of problem Equation (A8) by using any mathematical optimization method, like GA, SQP, etc., and compute the corresponding objective value y*; (e) Justify whether the design [x*, y*] satisfies the convergence condition Equation (A9). If not, add [x*, y*] into [X, Y], then recall steps (b)-(e), until the convergence condition is satisfied; If yes, optimization can finish: EI i y max − y min ≤ c (A9) From the procedure, it is said that solving Equation (A8) is not executed once, but executed again and again. Each solving needs to reconstruct the kriging model and EI function. Therefore, it reveals the nature that the MEI-based optimization method is a method of sequential modeling, indirectly optimizing the original problem.