Engineering with Computers (2006) 21: 289–295 DOI 10.1007/s00366-006-0018-x
ORIGINAL ARTICLE
Geng Tie Li Dequn Zhou Huamin
Three-dimensional finite element method for the filling simulation of injection molding
Received: Received: 28 December 2004 / Accepted: Accepted: 26 January 2006 / Published Published online: 20 May 2006 Springer-Verlag lag London London Limited 2006 Ó Springer-Ver
Abstract With the development of molding techniques, molded parts have more complex and larger geometry with nonuniform thickness. In this case, the velocity and the variation of parameters in the gapwise direction are considera considerable ble and cannot cannot be neglecte neglected. d. A three-dim three-dimenensional sional (3D) (3D) simula simulatio tion n model model can predic predictt the filling filling process more accurately than a 2.5D model based on the Hele–Shaw Hele–Shaw approximatio approximation. n. This paper gives gives a mathemathematical model and numeric method based on 3D model to perform more accurate simulations of a fully flow. The model model employs employs an equal-or equal-order der velocity– velocity–pres pressure sure interpolation method. The relation between velocity and pressure pressure is obtained obtained from the discretize discretized d momentum momentum equations in order to derive the pressure equation. A 3D control volume scheme is used to track the flow front. During calculating the temperature field, the influence of convection items in three directions is considered. The soft softwa ware re base based d on this this 3D mode modell can can calc calcul ulat atee the the press pressure ure field, field, veloc velocity ity field field and temper temperatu ature re field field in filling process. The validity of the model has been tested through the analysis of the flow in cavities. Keywords 3D Æ Equal-order interpolation Æ Injection molding Æ Simulation
1 Introductio Introduction n
During During injection injection molding, molding, the rheologic rheological al response response of polymer melts is generally non-Newtonian and nonisoG. Tie Æ L. Dequn Æ Z. Huamin State Key Laboratory of Mold & Die Technology, Huazhong Huazhong University University of Science and Technology, Technology, Wuhan, Hubei 430074, People’s Republic of China G. Tie (&) Machine Machine and Electric Electric Engineering Engineering College, Henan University of Technology, Technology, 450052 Zhengzhou, Henan, People’s Republic of China E-mail:
[email protected] Tel.: +86-0371-67758626 Fax: +86-372-3932808
thermal with the position of the moving flow front [1 [ 1– 3]. Because of these inherent factors, it is difficult to analyze the filling process. Therefore, simplifications are usually used. used. For example, example, in tradition traditional al middle-pl middle-plane ane model model and dual-domain model [4 [ 4, 5], the Hele–Shaw approximation [6 [6] is used. So both of these models are 2.5D models. In 2.5D model, the velocity and the variation of pressure in the gapwise direction are neglected except that the temperature is solved by FDM, and the filling of a mold cavity becomes a 2D problem in flow direction and a 1D problem in gapwise direction. As most of the inject injection ion molde molded d parts parts have have a sheetsheet-lik likee geomet geometry ry in whic which h the the thic thickn knes esss is much much smal smalle lerr than than the the othe otherr dimensions of the part, these models have been generally successful in predicting the advancement of melt fronts, pressure fields, and temperature distribution. The interest in 3D simulation of injection molding has increa increasedtreme sedtremendo ndous usly ly andsome progre progress ss has been been made made [7– 9] in the the past past few few year years. s. One One reas reason on is the the proc proces essi sing ng of larg largee and and comp comple lex x part parts. s. Wi With th the the deve develo lopm pmen entt of molding techniques, more and more molded parts have thick thick or nonun nonunifo iform rm thickn thickness ess,, such such as those those encoun encounter tered ed in gas gas-as -assis sisted ted inject injection ion moldin molding. g. In these these cases, cases, the velocity and the changes of parameters in the gapwise direction are considerable and cannot be neglected. On the other hand, the requirements on the performance of injection injection molded molded items items have been ever increasi increasing. ng. Several Several situations occurring during mold filling which cannot be accurate accurately ly predicted predicted using using the Hele–Shaw Hele–Shaw approxim approximatio ation n need to taken into account nowadays, such as the fluid behav behavior ior at the free free surfac surfacee (flow (flow front) front),, the fluid fluid behavi behavior or near and at the solid walls, the phenomenon occurring at merging merging of two or more fluid streams streams (weldlines) (weldlines),, and the kinematics in areas where shear and extensional deformations contribute significantly to the stress field (gates, ribs, etc.). A 3D simulation model should be able to generate complementary and more detailed information related lated to the flow flow charac character terist istics ics and stress stress distri distribut bution ionss in molded parts. This will be particularly important when dealin dealing g with with multic multicom ompon ponent ent mold mold filling filling and and with with molding of fiber-reinforced systems.
290
This paper presents a 3D finite element model to deal with the 3D flow of injection molding. In this model, the velocity in the gapwise direction is not neglected and the pressure also varies in this direction. An equal-order velocity–pressure formulation method [10– 12] is employed, and the relation between velocity and pressure is obtained from the discretized momentum equations. A 3D control volume scheme is introduced to track the flow front. During calculating temperature field, the influence of convection items in three directions is considered in order to get more exact results and to apply to the wider range of parts. Finally, the validity of the model has been tested through the analysis of some cases.
g0 ðT ; P Þ
.
1 þ ðg0 c_ sà Þ
1Àn
;
where n, c_ ; s* are non-Newtonian index, shear rate and material constant, respectively. Because there is no notable change in the scope of melt temperature during filling, Arrhenius model [13] for g0 is employed as following:
g0 ðT ; P Þ ¼ B exp
T b T
expðbP Þ;
where B, T b, b are material constants. 3 Finite element calculations for the pressure field
2 Governing equations
The pressure of melt is not very high during filling the cavity, so themelt is considered incompressible. Inertiaand gravitation are neglected as compared to the viscous force. With the above approximation, the governing equations, expressed in Cartesian coordinates, are as following: Momentum equations:
g¼
!
!
!
!
!
@ @ u @ @ v @ u @ @ w @ u 2g þ g þ þ g þ @ x @ x @ y @ x @ y @ z @ x @ z @ ðP Þ À ¼0 @ x @ @ v @ u @ @ v @ @ w @ v g þ þ þ g þ 2g @ x @ x @ y @ y @ y @ z @ y @ z : @ ðP Þ À ¼0 @ y @ @ w @ u @ @ v @ w @ @ w g þ þ g þ þ 2g @ x @ x @ z @ y @ z @ y @ z @ z @ ðP Þ À ¼0 @ z ð1Þ
!
3.1 Velocity–pressure relation In a 3D model, since the change of the physical quantities are not neglected in the gapwise direction, the momentum equations are much more complex than those in a 2.5D model. It is impossible to obtain the velocity–pressure relation by integrating the momentum equations in the gapwise direction, which is done in a 2.5D model. The momentum equations must be first discretized, and then the relation between velocity and pressure is derived from it. In this paper, the momentum equations are discretized using Galerkin’s method with bilinear velocity–pressure formulation. The element equations are assembled in the conventional manner to form the discretized global momentum equations and the velocity may be expressed as following: ui ¼ ~ ui À K iu
@ P @ P @ P ~ i À K iw vi ¼ ~ vi À K iv wi ¼ w ; @ x @ y @ z
where ~ ui ¼
X X X À
Axij uj À Bxij vj À C ijx wj =Axii ;
¼j i6
Continuity equation:
@ u @ v @ w þ þ ¼ 0: @ x @ y @ z
ð2Þ
Energy equation:
~ vi ¼
À
y
y
À
qC P
@ T @ T @ T @ T @ @ T K ¼ qC P u þv þw þ @ t @ x @ y @ z @ x @ x ; @ @ T @ @ T 2 þ þ ðK Þ þ gc_ K @ y @ y @ z @ z
y
C ijz wj À Az ij uj À Bz ij vj =C iiZ ;
¼j i6
y
Bij vj À Aij uj À C ij wj =Bii ;
¼j i6
~i ¼ w
the nodal pressure coefficients are defined as: K iu ¼
V
ð3Þ
where x, y, z are 3D coordinates and u, v, w are the velocity components in the x, y, z directions. P, T , q and g denote pressure, temperature, density and viscosity, respectively. Cross-viscosity model has been used for the simulations:
0Z @ 0Z @ 0Z @
K iv ¼
V
K iw ¼
V
1 A 1 A 1 A
N i dV =Axii ;
y
N i dV =Bii ;
N i dV =C iiz ;
ð4Þ
291
where Axij , Bxij , C xij , Ayij , Byj , C yij , Azij , Bzij , C zij represent global velocity coefficient matrices in the direction of x, y, z coordinate, respectively. K ui , K i , K w i denote the nodal pressure coefficients in the direction of x, y, z coordinate, respectively. The nodal values for K ui , K i , K w i are obtained by assembling the element-by-element contributions in the conventional manner. N i is element interpolation and i means global node number and j is, for a node, the amount of the nodes around it.
velocity field obtained by solving momentum equations does not satisfy continuity equation. The velocities are updated using the following relations:
3.2 Pressure equation
3.5 The tracing of the flow fronts
Substitution of the velocity expressions (4) into discretized continuity equation, which is discretized using Galerkin method, yields element equation for pressure:
The flow of fluid in the cavity is unsteady and the position of the flow fronts varies with time. Like in 2.5D model, in this paper, the control volume method is employed to trace the position of the flow fronts after the FAN (flow analysis network)[14]. But 3D control volume is a spacial volume and more complex than the 2D control volume. It is required that 3D control volumes of all nodes fill the part cavity without gap and hollow space. Two 3D control volumes are shown in Fig. 1.
v
v
Z
!
@ Ni @ Nk @ Ni @ Nk N j K ju P k þ N j K jv P k @ x @ x @ y @ y
V
@ Ni @ N k N j K jw P k þ dV @ z @ z @ Ni @ Ni @ Ni ~ j dV : ¼ N j ~ uj þ N j~ vj þ N j w @ x @ y @ z
Z
V
The element pressure equations are assembled in the conventional manner to form the global pressure equations.
3.3 Boundary conditions In the cavity wall, the no-slip boundary conditions are employed, e.g., u ¼ v ¼ w ¼ 0;
~ ~ ¼ 0; u¼~ v¼w
K iu ¼ K iv ¼ K iw ¼ 0
on an inlet boundary, u ¼ v ¼ w ¼ given
K iu ¼ K iv ¼ K iw ¼ 0:
3.4 Velocity update After the pressure field has been obtained, the velocity values are updated using new pressure field because the Fig. 1 3D control volumes. a Control volume of an internal node and b a boundary node
1 ui ¼ ~ ui À x Aii ~i À wi ¼ w
@ P N dV @ x
Z Z V
1 C iiz
N
1 vi ¼ ~ vi À y Bii
Z
N
@ P dV @ y
V
@ P dV : @ z
V
4 Finite element calculations for the temperature field
The temperature field plays an important role during injection molding process. Because the viscosity of the polymer varies with its temperature, so the variation of the temperature of polymer will have important influence to the injection molding process. Only after the temperature field during filling has been calculated exactly, the simulations for packing and cooling are meaningful. In 2.5D model, though the variation of the temperature in the gapwise direction is solved by FDM, the model is based on the Hele–Shaw approximation, which supposes the injection-molded parts are thin. As shown in Fig. 2, in 2.5D model, the triangular elements are meshed in the gapwise direction by creating finite difference grids, and the temperature in flow plane is represented by linear interpolation, and the temperatures in the gapwise direction are represented by FDM. In 2.5D models, the velocity in the gapwise direction is neglected, so only the thermal conduction item is considered in the gapwise direction. This paper gives a 3D
292
Thermal convection item and viscous heat item are anisotropic and has to do with the direction of flow. To keep the numerical stability, the upwind method is employed to handle the convection item and viscous heat item, e.g., only the contributions of the upriver elements from the nodes are considered when the convection item and viscous heat item are calculated. In the above equations, the time T is discretized using a forward-difference method: nþ1
@ T j T j ¼ @ t
Fig. 2 Illustrative finite difference in the gapwise direction
model for calculating the temperature field which considers the influence of convection items on three dimensions and suitable for the wider range of parts and has more exact results compared with the 2.5D models. According to the energy equation (3), by the use of Galerkin’s method, the equation for the temperature field can be expressed as following:
Z V
Z !
@ T N qC p dV ¼ @ t
N À qC p u
@ T @ T @ T þv þw @ x @ y @ z
V
@ @ T @ @ T @ @ T K K K þ þ þ þ gc_ 2 dV : @ x @ x @ y @ y @ z @ z Fig. 3 The test cavity. a The cavity dimension and b the meshed cavity
Fig. 4 Comparison between predicted shapes of flow front based on present 3D model ( a) and based on 2.5D model ( b). a Shape of 3D flow front and b shape of 2.5D flow front
À T jn
Dt
;
where Dt denotes time step. The element temperature equations are assembled in the conventional manner to form the global temperature equations. The overall procedure for pressure and temperature calculations is relaxation iterative. Because the pressure, velocity and temperature influence each other during the calculation, the temperature and pressure are coupled during the procedure.
5 Results and discussion
The first test cavity and dimensions are shown in Fig. 3a. The meshed 3D model of cavity is shown in Fig. 3b. The selected material is ABS780 from Kumbo. The parametric constants corresponding to the n,s*, B,T b and b of the five-constant cross-type viscosity model are 0.2638, 4.514 · 104 Pa, 3.13198043 · 10À7 Pa S, 1.12236 · 104 K, 0.000PaÀ 1. Injection temperature is 45°C, mold temperature is 250 °C, injection flow rate is 44.82 cu cm/s. ‘‘Fountain flow’’ is a typical flow phenomenon during filling. It has to do with the fluid near the center moving
293 Table 1 Material properties
Index
Material property
Unit
Reference value
1
Density (q) Specific heat (C p) Thermal conductivity (K ) Cross-type viscosity model N B T b b s
kg/m3
968.6
J/kg K
1.70
W/(m K)
0.140
Pa s K 1/Pa Pa
0.3783 1.0527 · 10À3 9.3841 · 103 0 1.955 · 103
2 3 4
Fig. 5 The example cavity
Fig. 6 The meshed cavity Fig. 7 The flow front at four different filling times. Time = 0.08, 0.36, 0.65 and 0.80 s
·
103
faster than the average across the thickness and upon catching up with the front, deflecting to move toward the walls, so the shape of the flow front is round like the fountain. In 2.5D models, the convection effects in the fountain region cannot be represented and the details of the fountain region are also lost, as shown in Fig. 4b. In presented 3D model, this fountain flow phenomenon can also be simulated. The round shape of the flow fronts at three filling times is illustrated clearly in Fig. 4a. Another example is typical of an industrial application as shown in Fig. 5. The outline dimensions of the cavity are 63.3 · 43.4 · 24.2 mm3 with a thickness 4 mm. ‘‘ ’’ represents the location of entrance.
294 Fig. 8 Temperature field on the plane Z = 13 in four filling times. Time = 0.08, 0.36, 0.64 and 0.82 s
The meshed cavity is as in Fig. 6. The injection temperature is 250°C, mold temperature is 45 °C, injection time is 0.82 s. The selected material is PS ASAHIPS 408. The material properties and the parametric constants corresponding to the five-constant cross-type viscosity model are specified in Table 1. Figure 7 shows the locations of flow fronts in four different filling times. A complex 3D flow field develops in the cavity and a rounded free surface is clearly seen. It can be seen that the filling process of the melt in the interior of the cavity can be predicted in the 3D model, and it is crucial for predicting more exactly the locations of the weldlines, the possibility of the air entrapments as well as the pressure and temperature distributions,
especially for these cavities with complex geometry and thick walls. Compared with the 2.5D model, which can only simulate the flow of melt in the surface of the cavity, 3D simulation model is suitable for the wider range of cavities and has more exact results. Figure 8 shows the temperature distributions on the plane Z = 13 in four filling times. It can be seen that there is a higher temperature in the interior of the cavity and the lower temperature near the cavity walls. The temperature near the entrance is even higher than the injection temperature due to the viscous heating of the melt: Figure 8 shows that a thermal layer is presented in the filled portion of the cavity, in which there is a variation of temperature from wall to the interior of cavity.
295
It can be seen that the heat transfer is mainly driven by convection and the conduction is rather small, so the thermal layer is very thin. It can be seen that the temperature distributions on the arbitrary section plane can be seen clearly in the 3D model. On the contrary, in the 2.5D model, only the average temperature in the thickness direction is shown on the middle-plane or the surfaces of the cavity, and for the thick or nonuniformthickness parts, which are not, suitable for the Hele– Shaw approximation, the results from 2.5D model have much error and even are mistakes.
6 Conclusion
A numerical model to simulate the filling of injection molding based on a 3D finite element model is presented in this paper. The 3D model uses the equal-order velocity–pressure formulation method and a 3D control volume scheme is adopted to track the flow front. During calculating temperature field, the influence of convection items in three directions is considered in order to get more exact results and to apply to the wider range of parts. Two parts have been employed as example to test the validity. It has been seen that 3D simulation model is suitable for the wider range of parts and has more exact results compared with 2.5D models. Acknowledgements The authors would like to acknowledge financial support from the National Natural Science Foundation Council of the People’s Republic of China, under Grant 20490220 and Research Foundation for PhD Candidates of Universities of the People’s Republic of China (20020487032).
References 1. Gotoh Terumasa, Iizuka, Miyamoto Masayuki et al (1986) Simulation of polymeric flows in the cavity filling process of injection molding. Sharp Tech J 34:63–70
2. Schlichting H (1968) Boundary-layer theory. McGraw-Hill, New York 3. Chen BS, Liu WH (1989) Numerical simulation and experimental investigation of injection mold filling with melt solidification. Polym Eng Sci 29:1039–1050 4. Dequn L (2002) New progress of flow simulation for plastic injection molding. China Int Forum Die Mould Technol 5:47– 48 5. Huamin Z, Dequn L (2002) Computer filling simulations of injection molding based on 3D surface model. Polym Plast Tech Eng 41:91–1021 6. Hiebert CA, Shen SF (1980) A finite-element/finite-difference simulation of injection molding filling process. J Non-Newtonian Fluid Mesh 7:1–32 7. Inoue Y, Higashi T, Matsuoka T (1996) Numerical simulation of 3-dimensional flow in injection molding. ANTEC Proc 1:744–748 8. Pichelin E, Coupez T (1998) Finite element solution of the 3D filling problem for viscous incompressible fluid. Comput Methods Appl Mech Eng 163:359–371 9. Hwang CJ, Kwon TH (2002) A full 3D finite element analysis of the powder injection molding filling process including slip phenomena. Polym Eng Sci 42:33–50 10. Prakash C, Patankar SV (1985) A control volume-based finiteelement method for solving the Navier–Stokes equations using equal-order velocity–pressure interpolation. Numer Heat Transfer 8:259–280 11. Prakash C (1986) An improved control volume finite-element method for heat and mass transfer, and for fluid flow using equal-order velocity–pressure interpolation. Numer Heat Transfer 9:253–276 12. Rice JG, Schnipe RJ (1986) An equal-order velocity–pressure formulation that does not exhibit spurious pressure modes. Comput Methods Appl Mech Eng 58:135–149 13. Hieber CA (1987) Injection and compression molding fundamentals. Marcel Dekker, New York 14. Tadmor Z, Broyer E, Gutfinger C (1974). Flow analysis method for solving flow problems in polymer processing. Polym Eng Sci 14:660–665