Journal of Materials Processing Technology 137 (2003) 74–77
Element modeling of FEM on the pressure field in the powder injection mold filling process Bing-yan Jiang*, Jue Zhong, Bai-yun Huang, Xuan-hui Qu, Yi-min Li Central South University, Changsha, Hunan 410083, P.R. China
Abstract In this paper, based on the procedure of derivation of the element equation, a simulation study is made in depth and the factors affecting the simulating precision are investigated. First the function of the fluidity from definition, physical significance and coupling solution aspects was analyzed in detail and a new expression, fluidity S, was presented. Then a viscosity model considering the effects of multi-components binder and powder loading was used. Finally, an algorithm for element kij of the stiffness matrix was set up, applying an element analysis method for non-Newtonian and non-isothermal feedstock melts in a thin mold cavity with symmetrical geometry. # 2002 Elsevier Science B.V. All rights reserved. Keywords: Powder injection molding; Pressure field; FEM; Element model
1. Introduction The frozen layer thickness is usually ignored in simulation of the injection mold filling processes. In this work, considering the frozen layer thickness, the physical model shown in Fig. 1 for non-Newtonian fluid under non-isothermal conditions is used. In this case, the simplified flowchart of the computer program for the filling process is shown in Fig. 2 [2]. There are four loops [3]: loop 1 is the pressure field iterative; loop 2 is the temperature field iterative; loop 3 is the flow front advancement at each time step until the cavity is filled; and loop 4 is a sub-loop, which determines the frozen layer thickness of the solving temperature. Obviously, the solution of the pressure field is the most internal loop: of course, it is most important. On the other hand, element analysis also is the core of the finite element analysis. Therefore, the element model of the pressure field was researched in depth in the following sections.
2. Pressure equation for the filling phase
2. The feedstock melt is assumed to be incompressible, i.e. the density r of the feedstock melt is constant. 3. The thermal conductivity k of the feedstock melt is assumed to be constant. 4. The flow field is assumed to be symmetrical about the center line of the cavity. 5. Gravity and the force of inertia of the feedstock melt are neglected. Using an assumption about the feedstock melt behavior and the geometric symmetry of the thin cavity, the following governing equations for the liquid phase [2,5] can be obtained: @u @v @w þ þ ¼0 @x @y @z @p @ @u Momentum equation : ¼ Z ; @x @x @z @p @ @v @p ¼ Z ¼0 ; @y @y @z @z @T @T @T Energy equation : rcp þ vx þ vy @t @x @x
Continuity equation :
Basic assumptions:
¼ Z_g2 þ k
1. The feedstock melt is assumed to be a continuum.
* Corresponding author. E-mail address:
[email protected] (B.Y. Jiang).
where 2 2 @vx @vy 2 þ g_ ¼ @z @z
0924-0136/02/$ – see front matter # 2002 Elsevier Science B.V. All rights reserved. PII: S 0 9 2 4 - 0 1 3 6 ( 0 2 ) 0 1 0 7 0 - 1
@2T @z2
(1)
(2)
(3)
B.Y. Jiang et al. / Journal of Materials Processing Technology 137 (2003) 74–77
75
where the viscosity is not too high and z is not too small. Of course, this region is determined by the thickness of the frozen layer. It follows that the calculation of the fluidity and that of the frozen layer thickness are related. Clearly S differs at different z, but the distribution of S along the z-direction can be obtained. A new more appropriate definition of fluidity is R hz presented by sðzÞ ¼ 0 ð~z2 =ZÞ d~z. Physically, the fluidity determines the ease with which the melt can be forced through the mold cavity. The fluidity gives a measure of the ease with which the feedstock melt flows. A high value of fluidity is associated with easy flowing materials. Conversely, a low value of fluidity indicates that the feedstock melt offers a high resistance to flow [2]. From the aspect of the solution of the governing equation, fluidity S is a key parameter to solve the pressure equation. When a viscosity model related to temperature and shear rate is used, the pressure algebraic equation discretized is linear. When a model related to pressure is used, the pressure algebraic equations are non-linear, while the temperature field is calculated through the energy equation, S is a crucial parameter coupling the pressure equation with the energy equation. Studies on the evaluation of fluidity are very rare, thus this topic needs to be studied further.
Fig. 1. The frozen layer thickness.
4. The constitutive equation of the feedstock
Fig. 2. Simplified flowchart of the program.
Further, the energy equation for the solid phase is rs cps
@Ts @ 2 Ts ¼k 2 @t @zs
(4)
It is important for the study on feedstock rheology to improve the precision of simulation in powder injection molding. According to the literature, there are over 100 viscosity models. The selection of a practical viscosity model for simulation is much more important. Here the viscosity model was used in the following form [1]: n X Binder viscosity : ln Zb ¼ wi ln Zi (7) i¼1
f 2 Zf ¼ Zb 1 fm E 1 1 Pure binder : Z ¼ Z0 exp k T T0
Integrating the continuity equation and the momentum equation, the following equation is obtained: Z h 2 @ @p @ @p z S S dz (5) þ ¼ 0; where S ¼ @x @x @x @y 0 Z
Feedstock viscosity :
Applying a two-dimensional divergence and gradient operator to Eq. (5):
where Zi indicates the viscosity of component i (polymer), wi the weight of the component i (polymer), Zb the viscosity of the binder, Zf the viscosity of the feedstock, f the powder loading of the feedstock, fm the powder maximum loading of the feedstock, E the activation energy for viscous flow, K the Boltzmann’s constant, and Z and Z0 are the viscosities of the pure binder at absolute temperature T and reference temperature T0, respectively.
divðSrpÞ r Srp ¼ 0
(6)
3. Coupling parameter—fluidity, S From Eq. (5), it is seen that the greatest contribution to S comes from the region where the viscosity starts to increase rapidly. This is because around the center line, z is small and so the contribution to S is small, whereas at the wall, where z is large, the viscosity is extremely large and again the contribution to S is small. The most significant region is therefore
(8)
(9)
5. Element model on pressure equation An algebraic equation using node pressure as the unknown variable was set up, as for any application problem
76
B.Y. Jiang et al. / Journal of Materials Processing Technology 137 (2003) 74–77
with FEM, to obtain the following equation: ½K e fPge ¼ fQge
(10)
where [K]e is called the element stiffness matrix, {P}e the nodal pressure vector, {Q}e the nodal flow rate vector, and the cavity mold is meshed with triangular elements. The above equation can be expressed in detail as 2 38 9 8 9 k11 k12 k13 > = > = < p1 > < Q1 > 6 7 ¼ (11) k k k p Q2 4 21 22 23 5 2 > ; > ; : > : > k31 k32 k33 p3 Q3 All the elements have the same form, there being only different element values in the element stiffness matrix. Hence one of the core factors in finite element analysis is element analysis. The element model is essentially the determination of element arithmetic at the element stiffness matrix. The assembly of the system stiffness matrix is another problem, related to the pressure and temperature boundary conditions at the solution domain. The key point of the former is element arithmetic at the stiffness matrix [K]e, whilst the key point of the latter is the application of the boundary conditions and the assembly of the system matrix. Although these two parts are indiscerptible, here at first the element modeling of FEM on the pressure field for the powder injection mold filling process is discussed. A linear interpolation function is used for the pressure within an element [5]. Hence pressure Pðx; yÞ at any point within an element is approximately expressed by the shape function or area coordination [2–4]: ~ pðeÞ ðx; yÞ ¼
3 X
Ni pi ¼
i¼1
3 X
xi p i
(12)
i¼1
where Ni is a linear shape function, pi the pressure at node i, and xi is the area coordinate on a triangular element. It can be shown that shape function Ni in the linear interpolation function is the area coordinates N i ¼ xi
(13)
where S is the boundary of the element. Rearranging this equation Z Z Sr~p rNi dA ¼ Ni Sr~p ~ n dS ði ¼ 1; 2; 3Þ (16) S
A
r~p; rNi were expressed in area coordinates and substituted into Eq. (16) to enable three algebra equations to be obtained. When i ¼ 1, the LHS of Eq. (16) is Z Z S fðy223 þ x232 Þp1 Sr~p rN1 dA ¼ S r~p rN1 dA ¼ 4A A A þ ðy31 y23 þ x13 x32 Þp2 þ ðy12 y23 þ x21 x32 Þp3 g
The solution of the right-hand side (RHS) of Eq. (16) is the line integration of the triangular element boundary, in Fig. 3, ~ n is the unit outward normal to the element side, and the element boundary S is composed of the three lines S1 ; S2 and S3 . r~p ~ n is actually the direction derivative @p=@n of the function P on the element boundary S. Then the flow rate of the nodes and the flow rate in unit length along the side can be obtained: Q1 ¼ 12 q2 jS2 j þ 12 q3 jS3 j;
Q2 ¼ 12 q3 jS3 j þ 12 q1 jS1 j;
Q3 ¼ 12 q2 jS2 j þ 12 q1 jS1 j q1 ¼ S
@p ; @n1
q2 ¼ S
(18) @p ; @n2
q3 ¼ S
@p @n3
(20)
A
(14) From the divergence equation of the function product, Eq. (14) can be expressed as Z Z r ðNi Sr~ pÞ dA Sr~ p rNi dA ¼ 0 (15) A
A
Applying the Gauss theorem in the plane to the first term on the left-hand side (LHS) of Eq. (15): Z Z Ni Sr~ p ~ n dS Sr~ p rNi dA ¼ 0 ði ¼ 1; 2; 3Þ S
A
(19)
where jS1 j; jS2 j; jS3 j indicate the length of the triangle element boundary S1, S2, S3, respectively, and q1, q2, q3 are the unit length flow rate along the side S1, S2, S3, respectively. Q1, Q2, Q3 are the volume flow rates at nodes 1, 2, 3, respectively. For the RHS of Eq. (16), Z Z Z @p @p dS ¼ S N1 Sr~p ~ n dS ¼ x1 S x dS @n @n2 S2 1 S S2 [S3 Z @p 1 1 x1 dS ¼ q2 jS2 j þ q3 jS3 j þS @n3 S3 2 2
Applying the Galerkin method to Eq. (6): Z Z Ni Rð~ pÞ dA ¼ Ni ðr Sr~ pÞ dA ¼ 0 ði ¼ 1; 2; 3Þ A
(17)
Fig. 3. Parameters of the element.
B.Y. Jiang et al. / Journal of Materials Processing Technology 137 (2003) 74–77
where substituting for the LHS and RHS of Eq. (16): S fðy223 þ x232 Þp1 þ ðy31 y23 þ x13 x32 Þp2 4A þ ðy12 y23 þ x21 x32 Þp3 g ¼ Q1
where xij ¼ xi xj ;
yij ¼ yi yj ;
(21)
(22)
S fðy23 y12 þ x32 x21 Þp1 þ ðy31 y12 þ x13 x21 Þp2 4A þ ðy212 þ x221 Þp3 g ¼ Q3
(23)
The 3 3 matrix on the left is called the element stiffness matrix. It is symmetric, therefore to describe the element stiffness matrix, only six terms are required, in which S 2 S ðy þ x232 Þ; ðy31 y23 þ x13 x32 Þ; k12 ¼ 4A 23 4A S S 2 ðy23 y12 þ x32 x21 Þ; ðy þ x213 Þ; ¼ k22 ¼ 4A 4A 31 S 2 S ðy þ x221 Þ; ðy12 y31 þ x21 x13 Þ (25) ¼ k23 ¼ 4A 12 4A
k11 ¼
k33
Considering symmetry : kij ¼ kji
x1 x2 x3
y1 y2 y3
6. Conclusions
The above set of three equations (21)–(23) are called the element equations. They give the relationship between pressure p1, p2, p3 and flow rate Q1, Q2, Q3 for the triangular element written in matrix notation as follows: 2 3 y223 þ x232 y23 y31 þ x32 x13 y23 y12 þ x32 x21 S 6 7 y231 þ x213 y12 y31 þ x21 x13 5 4 y23 y31 þ x32 x13 4A y23 y12 þ x32 x21 y12 y31 þ x21 x13 y212 þ x221 8 9 8 9 > = > = < p1 > < Q1 > (24) p2 ¼ Q2 > ; > ; : > : > p3 Q3
k13
1 1 A ¼ 1 2 1
(27)
In Eq. (16) let i ¼ 2; 3, applying the same way, to obtain S fðy23 y31 þ x32 x13 Þp1 þ ðy231 þ x213 Þp2 4A þ ðy12 y31 þ x21 x13 Þp3 g ¼ Q2
77
(26)
1. A more appropriateR definition of fluidity S was hz presented, i.e. sðzÞ ¼ 0 ð~z2 =ZÞ d~z. 2. The coupling of the pressure equation with the energy equation is decoupled by a crucial parameter, fluidity S. Whether the pressure algebraic equations are linear or non-linear is determined by the relation between S and P. 3. An algorithm for the element kij of stiffness matrix is derived. 4. The analysis process on the element model paid the foundation for abstracting nodes class and element class in Cþþ in an object oriented analysis. Acknowledgements This work was supported by grants from the Lab Open Fund of The National Key Laboratory for Powder Metallurgy, partially from 973 projects of G2000067200. References [1] Y. Li, B. Huang, X. Qu, Viscosity and melt rheology of metal injection moulding feedstock, Powder Metall. 42 (1) (1999) 86–90. [2] P. Kenndey, Flow Analysis of Injection Molds, Hanser Publishers, Munich, 1995. [3] C.A. Hieber, S.F.A. Shen, Finite element/finite difference simulation of the injection molding filling process, J. Non-Newtonian Fluid Mech. 7 (1980) 1–8. [4] L.A. Najmi, D. Lee, Modeling of mold filling process for powder injection molding, Polym. Eng. Sci. 31 (15) (1991) 1137–1139. [5] B.S. Chen, W.H. Liu, Numerical simulation and experimental investigation of injection mold filling with melt solidification, Polym. Eng. Sci. 29 (15) (1989) 1041–1043.