A new explicit simulation for injection molding and its validation.

INTRODUCTION

Injection molding is a widely used technology in manufacture. For the process simulation of polymer injection molding of limited thickness models (1), (2), Hele-Shaw model is generally adopted. The tracking of filling front is usually modeled by a pseudo-concentration method (3). However, 3D model is furthermore required when the model becomes more complicated than a quasi 2D one (4), (5). Moreover, for the simulation of bi-injection, coinjection, and gas-assisted injection processes, 3D models represent the obvious advantages (6). Adequate methods of simulation combined with experiments are the valid ways to predict the product quality and determine the process parameters in injection molding (7-9). To improve the efficiency of computation for mold filling problem, many efforts have been made. The performance to solve the industrial problems remains still a challenge when complicated models are involved. The explicit algorithms proposed in the literature for mold filling (10), (11) during the last 10 years make a significant advance by virtue of the application of dynamic models and the associated lumped mass matrices. In these explicit simulations, MINI elements are traditionally employed. These simulations are much faster than the implicit ones and their validity is proven by the experiments. Nevertheless, global solutions are always necessary to evaluate the pressure fields at each time step. It results in the increase of computational cost exponentially proportional to the DOF number. It is a serious obstacle for industrial applications, in which the 3D large-scale problems are always involved.

On the basis of previous works (12), (13), a new explicit algorithm is proposed and tested for simulation of the mold filling. The most important feature is its purely vectorial strategy. Except for the operations at element level, neither global solution nor the construction of global matrix is required. This method provides a powerful approach to predict the filling process. The computational cost is about linearly proportional to the DOF number. The incompressibility is maintained by a special feedback regulation strategy. Moreover, the proposed algorithm is very easy to be parallelized for high performance computation on multiclusters, as there are only vectorial operations performed in uncoupled manner. The comparison of the present numerical results with previous simulations proves the validity and efficiency of the newly proposed algorithm.

MECHANICAL MODELING OF MOLD FILLING

Because of the nature of filling process, Eulerian description is used in modeling of the mold filling. Let t be the current instant in injection process: t [member of] [0, [t.sub.fn], [t.sub.fn] represents the final time corresponding to the end of filling process. The whole position in the mold, each represented by the spatial vector X, is defined as a set [OMEGA]. This set consists of two different parts at each instant: the part filled by polymer and the remained void part, practically filled by air. They are represented by two adjoined subsets [[OMEGA].sup.F] (t) and [[OMEGA].sup.V](t). A field variable F(X,t) is used to describe the filling state, which takes value I in the filled portion and value 0 in the void portion.

The injection flow in filled part of the mold cavity is supposed to be incompressible. It should obey the momentum and mass conservation equations. When Reynolds number is small, as the case of mold filling, the momentum conservation is represented by a Stokes equation. This equation is expressed, respectively, in the filled and void portion as:

[for all] X [member of][[OMEGA].sup.F], [[rho].sub.f] [[partial derivative]V]/[[partial derivative]t] = - [nabla]P + [nabla]*([[mu].sub.f][[epsilon].sub.f]) + [[rho].sub.f]g (la)

[for all] X [member of] [[OMEGA].sup.v], [[rho].sub.a][[partial derivative]V]/[[partial derivative]t] = - [nabla]P + [nabla]*([[mu].sub.a][[epsilon].sub.a]) + [[rho].sub.a]g (lb)

where V is the velocity vector: P represents the pressure field; [[mu].sub.f] and [[mu].sub.a] the viscosities in two different mold parts; [[rho].sub.f] and [[rho].sub.a] their apparent densities in different parts; [[epsilon].sub.f] and [[epsilon].sub.a] are the shear rate in filled and void part; and g the gravity acceleration. To keep the numerical stability, Active density is applied for [[rho].sub.a] in the void mold part.

The constitutive relations for filled and void portion take generally the following form, respectively.

[[sigma]'.sub.f] = 2[[micro].sub.f](T, [[epsilon].sub.f])[[epsilon].sub.f] (2a)

[[sigma]'.sub.a] = 2[[micro].sub.a](T, [[epsilon].sub.a])[[epsilon].sub.a] (2b)

where [[sigma]'.sub.f] and [[sigma]'.sub.a] are the deviatoric Cauchy stress tensors in filled and void portion; T indicates the temperature value; and [[epsilon].sub.f] and [[epsilon].sub.a] are the equivalent values of shear rate deviators.

The constitutive relation is strongly dependent on the local temperature values, so that the variation of temperature field in injection molding is very important to determine the variation of viscosity. Then the evaluation of temperature field during the simulation of mold filling is necessary. Once the polymer is considered as isotropic and the Fourier model for heat flux is applied, the heat transfer is controlled by the following advective-diffusive equation:

[for all] X [member of] [[OMEGA].sup.F], [[rho].sub.f][C.sub.f]([[partial derivative]T]/[[partial derivative]t] + V * [nabla]T) = [nabla] * ([k.sub.f][nabla]T) + [[sigma]'.sub.f]: [[epsilon].sub.f] (3a)

[for all] X [member of] [[OMEGA].sup.v], [[rho].sub.a][C.sub.a]([[partial derivative]T]/[[partial derivative]t] + V * [nabla]T) = [nabla]*([k.sub.a][nabla]T) + [[sigma]'.sub.a]:[[epsilon].sub.a] (3b)

where [C.sub.f] and [k.sub.f] are, respectively, the specific heat and thermal conductivity coefficient of filled polymer, [C.sub.a] and [k.sub.a] are, respectively, the specific heat and thermal conductivity coefficient of air, [[sigma]'.sub.f]:[[epsilon].sub.p] and [[sigma]'.sub.a]:[[epsilon].sub.a] stand for the dissipations associated to viscous polymer flow and viscous air flow.

It should be mentioned that the incompressibility condition is attentively verified only in the filled part, although the same operation is performed in the void part to keep simplicity of the algorithm. By mass conservation, one has the following equation:

[for all] X [member of] [[OMEGA].sup.F], [nabla] * V = 0 (4)

At each instant t during the injection course, the evolution of filling state variable is governed by an advection equation:

[for all] X [member of] [OMEGA], [[partial derivative]F]/[[partial derivative]t] + V * [nabla]F = 0 (5)

SOLUTION PROCEDURE WITH FULLY VECTORIZED EXPLICIT STRATEGY

Instead of the use of MINI elements or reduced integration (14), the proposed algorithm is based on the elements with equal order interpolation. A smoothing procedure is performed systematically without piecewise operations. So the interpolations of different field variables are simply written as:

{V, P, F, T} = [m.summation over (i = 1)][N.sub.i]{[V.sub.i], [P.sub.i], [F.sub.i], [T.sub.i]} (6)

in which [N.sub.i] stands for the interpolation functions and m is the node number in an element.

In the explicit algorithm, the momentum conservation Eqs, la and lb; can be spitted into two fractional steps: (1) solution for the viscous diffusion and (2) the maintaining of incompressibility conditions.

Viscous Diffusion

The fractional step for viscous diffusion is expressed as:

[for all] X [member of] [[OMEGA].sup.F], [[rho].sub.f] = [[partial derivative]V]/[[partial derivative]t] = [nabla] * ([[mu].sub.f][[epsilon].sub.f]) + [f.sup.ext] (7a)

[for all] X [member of] [[OMEGA].sup.v], [[rho].sub.a] = [[partial derivative]V]/[[partial derivative]t] = [nabla] * ([[mu].sub.a][[epsilon].sub.a]) + [f.sup.ext] (7b)

where [f.sup.ext] stands for the external loads.

Equations 7a and 7b are discretized by the standard GFEM procedure, it results in:

M [([V.sub.n+] - [V.sub.n])]/[[DELTA]t] = [F.sup.[sigma]] ([V.sub.n]) + [F.sup.ext] (8)

In above equations, [V.sub.n+] represents the intermediate values for solution of the velocity fields, [V.sub.n] is the velocity at [t.sub.n] instant, [F.sup.[sigma]] represents the effects of viscous diffusion. [F.sup.ext] is the external load vector associated to gravity contribution and boundary conditions, and M is a lumped mass matrix built as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (9)

in which A stands for an operator for FEM assemblage and N is the matrix of interpolation function at element level. The mass matrix is lumped into diagonal form to get the solution uncoupled and explicit.

Incompressibility

The fractional step for incompressibility is expressed as:

[for all] X [epsilon] [[OMEGA].sup.F], [[rho].sub.f] [[partial derivative]V]/[[partial derivative]t] = -[nabla]P (10a)

[for all] X [epsilon] [[OMEGA].sup.v], [[rho].sub.a] [[partial derivative]V]/[[partial derivative]t] = -[nabla]P (10b)

It is mentioned that the solution of Eq. 10b is only for maintaining integrity of the model solution, rather than keeping incompressibility in the void parts.

In the previous algorithm using MINI element, Eqs. 10a and 10b are discretized by the standard GFEM procedure, which results in Eq. 11 to solve the pressure held. With the obtained pressure field, the intermediate velocity field is corrected by Eq. 12 to meet incompressibility.

[for all] X [member of] [OMEGA], ([[~.G].sup.T][M.sup.-1][~.G])P = 1/[DELTA]t [[~.G].sup.T][V.sub.n+] (11)

[for all] X [member of] [OMEGA], M ([V.sub.n + 1] - [V.sub.n+])/[DELTA]t = - [~.G]P (12)

where M is the mass matrix with the use of MINI elements, [V.sub.n + 1] is the velocity at [t.sub.n+1] instant, [~.G] is the gradient operator, and [[~.G].sup.T] is the divergent one defined in the global manner. They are constructed in the following forms:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (13)

At each Gauss point, [D.sup.+] is defined as a divergence operator [nabla] * V = [D.sup.+] [V.sup.e], which is constructed on the basis of interpolation functions. The term [[~.G].sup.T][M.sup.-1] [~.G] in Eq. 11 is a matrix, and so Eq. 11 should be solved globally. If it is supposed that the void portion possesses the same mass density as in the filled one, the term [[~.G].sup.T][M.sup.-1][~.G] becomes a constant matrix. This matrix is decomposed immediately in Choleski form after construction and used repeatedly in all the simulation time steps. However, the global solutions of pressure fields in each time steps are the most expensive ones in the previous explicit algorithms. Thousands of the time steps result in the tremendous number of global solutions for the pressure fields, which conducts the simulation of large scale problems to be extremely heavy charges in computation. In fact, people need to know only the pressure fields at some significant prescribed instants. The output of pressure fields for all the steps is not necessary.

To remedy such a drawback, a new method is proposed for the satisfaction of incompressibility. The intermediate velocity field obtained by the fractional step for the viscous diffusion effect is then corrected by a special factor [lambda]. This factor is determined by the respect of incompressibility. It is a local operation in which each localized correction at element is proportional to the deviated value regarding to incompressibility in the field of intermediate velocity. However, a smoothing operation is necessary for the deviated values to avoid numerical instabilities, as only the elements of equal interpolation order are used. The most important advantage of the new vectorial algorithm remains in the fact that the solution of pressure fields is no more necessary during the simulation. However, the pressure fields can be achieved by postprocessing of the new algorithm. One can prescribe some significant instants in the simulation to stock the concerned data, which includes the discretized values of pressure's gradient fields. Incorporating with the boundary condition on pressure, the pressure field can be solved globally by a special module of the postprocessing.

To satisfy the incompressibility condition, the following equations must be dealt to correct the intermediate velocity fields:

[for all] X [member of][[OMEGA].sup.F], [[rho].sub.f] [[partial derivative]([V.sub.n + 1] - [V.sub.n+])]/[[partial derivative]t] = - [nabla][[lambda].sup.f] [([nabla].[V.sub.n+])I] (14a)

[for all] X [member of] [[OMEGA].sup.v], [[rho].sub.a] [[partial derivative] ([V.sub.n + 1] - [V.sub.n+])]/[[partial derivative]t] = - [nabla][[lambda].sup.a] [([nabla].[V.sub.n+])I] (14b)

In above equations, incompressibility condition in the filled part is accounted by the right hand term with scale factor [[lambda].sup.f].

These Eqs. 14a and 14b are discretized by GFEM. The incompressible velocity field [V.sub.n + 1] is obtained directly from:

M ([V.sub.n + 1] - [V.sub.n+])/[DELTA]t = [F.sup.p]([V.sub.n+]) (15)

where [F.sup.p] stands for the term of correction to satisfy the incompressibility condition. This term is built as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (16)

where [[bar.d].sub.sm]([V.sub.n+]) is the deviated values at each Gauss points caused by intermediate velocity field [V.sub.n+], which is smoothened locally in the patch with its neighbor elements, B is the derivative matrix of the interpolation function N used to calculate the deformation rate in each element. [I.sup.c] is a constant matrix, and [[lambda].sub.n + 1] is the factor of regularization for filled and void domains, respectively, in the mold.

The deviated value d for the intermediate field [V.sub.n+] at each Gauss point is defined as:

d([V.sub.n+]) = 1/3 tr([[epsilon].sub.n+]) (17)

where [[epsilon].sub.n+] is the shear rate calculated by intermediate velocity field [V.sub.n+].

In mixed formulations using MINI element for incompressible flows, the interpolation of the velocity field is one order higher than the pressure field. Instead of it, the reduced integration can be used. As an alternative method, a smoothing operation is chosen in the new proposed algorithm to avoid the numerical instabilities. It is a systematic operation without piecemeal treatments (15). The nodal values of the deviated term [d.sub.sm] are obtained by the relationship:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (18)

in which [M.sub.0] is a pseudo lumped mass matrix in diagonal form, built as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (19)

So Eq. 16 can be written as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (20)

with the relationship [[bar.d].sub.sm] = N[d.sud.sm.sup.e] at Gauss points.

The same operations are necessary in the void mold part to maintain integrity of the model solution. As it is not the part in which we are really interested, one can use fictive factors to avoid the numerical instabilities.

Determination of the Parameter [lambda].

The determination of coefficient [[lambda].sub.n + 1] is a key procedure in new method. The average deviated value in the filled part [[OMEGA].sub.F] is used to determine the [[lambda].sub.n + 1] value for filled part in each time step. The way to determine [[lambda].sub.n+1] is expressed as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (21)

where the notation [<>.sub.sum] denotes the sum of all nodal values in the filled part, and D is a divergence operator. When [[lambda].sub.n + 1] is obtained, it can be introduced in Eqs. 15 and 16 to verify the incompressibility condition.

Determination of Filling State

Equation 5 is discretized by the standard Taylor-Galerkin procedure, it results in:

[M.sub.0][([F.sub.n + 1] - [F.sub.n])]/[[DELTA]t] = [[K.sup.a] + [K.sup.ad]][F.sub.n] (22)

where [F.sub.n] and [F.sub.n + 1] are the filling state vectors at [t.sub.n] and [t.sub.n + 1] instant, respectively. [K.sup.a] presents advection effect, and [K.sup.ad] presents diffusion effect.

Overcoming the Instability Caused by Discontinuous Mesh Filling

Under the frame of Eulerian description, the filling state is defined element by element for the integration of different behaviors. The filling of finite element model is in fact discontinuous because each element should be judged into two extreme states: fully filled or fully empty. This discontinuity is generally not significantly disturbing because the element sizes are generally enough fine for the simulation by an explicit algorithm.

However, for the application of the new algorithm, such an effect may induce the severe instability at the beginning of filling simulation. As an algorithm specially designed for simulation of the filling process, the incompressibility is retained only in filled portion with Eulerian mesh. The incompressible flow in filled domain is updated at each time step on the basis of previous time step. Then the filling should advance gradually to let the algorithm running.

While in the starting stage of filling simulation, the filling ratio changes suddenly compared with their previous states. As only one or few layers of the elements are filled at the beginning of simulation, the filling of one more-layer will cause significant discontinuity for the update of filled domain. Such an important discontinuity will cause instability for the fractional step to keep incompressibility in the filled domain, as the filling is not in a gradual manner.

To overcome this instability, a smoothing step is necessary for the determination of dilatancy parameter [lambda] at the starting stage of simulation. This parameter should be limited to prevent the changes in a violent manner because of the sharp change of filled portion. The following equation can be used to remove such a type of instability:

[[lambda].sub.n + 1] = [[lambda].sub.n] + [R.sub.f][[epsilon].sub.s] ([[lambda]'.sub.n + 1] - [[lambda].sub.n]) (23)

in which [[lambda].sub.n + 1] is the dilatation parameter used for the present time step and [[lambda].sub.n] is the one for last time step. [[lambda]'.sub.n + 1] is the parameter determined by feedback procedure, as expressed in Eq. 21. [R.sub.f] stands for the filling ratio with respect to total volume of the meshed model. [R.sub.f][member of] (0,1), [R.sub.f] = 1 represents the fully filled state. Generally, this factor takes a very small value at the beginning of simulation, as often one layer of the elements on inlet is supposed to be filled for the first time step. [[epsilon].sub.s] is a smoothing coefficient that takes a value significantly less than one.

This smoothing procedure is very effective to keep stability of the algorithm against the effect of discontinuous filling at the starting stage of simulation. As the time steps are chosen sufficiently small to fill one layer elements in several time steps, the incompressibility can be well kept in the filled domain by the stable corrections on intermediate velocity fields.

The Calculation of Heat Transfer and Temperature Field

The heat transfer in Eqs. 3a and 3b is dealt with both advection and diffusion effects. These two effects are solved usually by an implicit solution but subjected to the limit on time steps (16). A fractional step method proposed by Lewis et al. (17) permits to evaluate the temperature field explicitly. This method is retained in the present work as it provides a very good efficiency. The fractional advection and diffusion steps are written as:

[for all] X [member of] [OMEGA], [[partial derivative]([T.sub.n+] - [T.sub.n]]/[[partial derivative]t] = - V * [nabla]T (24)

[for all] X [member of] [[OMEGA].sup.F], [([T.sub.n + 1] - [T.sub.n+])]/[[partial derivative]t] = 1/[[rho].sub.f][C.sub.f] ([nabla] * ([k.sub.f][nabla][T.sub.n+]) + [[sigma]'.sub.f]:[[epsilon].sub.f]) (25a)

[for all] X [member of][[OMEGA].sup.v], [([T.sub.n + 1] - [T.sub.n+])]/[[partial derivative]t] = 1/[[rho].sub.a][C.sub.a] ([nabla] * ([k.sub.a][nabla][T.sub.n+]) + [[sigma]'.sub.a]:[[epsilon].sub.a]) (25b)

where [T.sub.n] and [T.sub.n + 1] are the temperature fields at instant [t.sub.n] and [t.sub.n + 1], and [T.sub.n + 1] is an intermediate temperature field used in the fractional steps method.

Same as for the evaluation of filling state by Taylor-Galerkin method, Eq. 24 can be discretized in following form:

[for all] X [member of][OMEGA], [M.sub.0] [([T.sub.n+] - [T.sub.n])]/[[DELTA]t] = - [[K.sup.ad](V) + [K.sup.df](V)][T.sub.n] (26)

By Galerkin method, the solution of Eqs. 25a and 25b can be written in the following manner:

[for all] X [member of] [OMEGA], [M.sub.0] [([T.sub.n + 1] - [T.sub.n+])]/[[DELTA]t] = - [[K.sup.th][T.sub.n+] + Q + [q.sup.c] (27)

where Kth is the thermal diffusion term, Q is the heat dissipation term, and [q.sup.c] is the heat convection term.

Fractional Steps in Both Algorithms

To compare the solution strategies of Stokes equation between the new vectorial algorithm and the previous one using MINI element, the fractional steps of both algorithms are shown in Table 1. One can see that the important difference between both algorithms locates on the fractional step to satisfy incompressibility condition. In previous algorithm, a global solution scheme must be performed to get the pressure fields. The global solution costs a large amount of computation when large model is involved. The cost is about exponentially proportional to the degree of freedom number. This is a serious obstacle for practical industrial applications.

TABLE 1. Comparison of the fractional steps in two algorithms.

1. Viscous diffusion      [for all] [CHI] [member of] [OMEGA] M
(local solution)          [(V *- [V.sub.n])]/[[DELTA]t] =
                          [F.sup.[sigma]]([V.sub.n]) + [F.sup.ext]

2(a). Incompressibility   (1) Coefficient [[lambda].sup.F] (local
condition realized in     solution)
the new algorithm (local
solution)

                          [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE
                          IN ASCII]

                          (2) Velocity fields (local solution)

                          [for all] [CHI] [member of] [OMEGA] i
                          [member of] {F, V}, M [([V.sub.n+1] - V *]/
                          [[DELTA]t]

2(b). Incompressibility   (1) Pressure held (global solution)
condition realized in
the previous algorithm
(including global
solution)

                          [for all] [CHI] [member of] [OMEGA],
                          ([[~.G].sup.T] M [.sup.I][~.G])P =
                          1/[DELTA]t [[~.G].sup.T] V *

                          (2) Velocity fields (local solution)

                          [for all] [CHI] [member of][OMEGA], M
                          [([V.sub.n+1] - V *)]/[[DELTA]t] = [~.-G]P

3. Filling state (local   [for all] [CHI] [member of] [OMEGA],
solution)                 [M.sub.0] [([F.sub.n+1] -
                          [F.sub.n])]/[[DELTA]t] = [[K.sup.a] +
                          [K.sup.ad]][F.sub.n]

4. Heat transfer (local   (1) Advection
solution)

                          [for all] [CHI] [member of]
                          [OMEGA], [M.sub.0] [([T.sub.n+] -
                          [T.sub.n])]/[[DELTA]t] = -[[K.sup.ad](V) +
                          [K.sup.ad] (V) + [K.sub.df] (V)]
                          [T.sub.n]

                          (2) Viscous diffusion
                          [for all] [CHI]
                          [member of] [OMEGA], [M.sub.0]
                          [([T.sub.n+1] -[T.sub.n+])]/[[DELTA]t] =
                          [K.sup.th][T.sub.n+] + Q + [q.sup.c]

Note: (a) new vectorial algorithm and (b) previous explicit algorithm.

The new algorithm is designed especially for the simulation of mold filling process. The incompressibility condition is maintained by feedback corrective scheme. This method is totally vectorized in numerical operation except for the small matrices at element level. So the computational cost is largely reduced. Moreover, the proposed algorithm is very easy to be parallelized for high performance computation on multiclusters, as there are only vectorial operations performed in uncoupled manner. This work builds also a solid basis for the development of biphasic vectorial algorithm, which can predict the important segregation effect efficiently in metal injection molding.

NUMERICAL EXAMPLES AND COMPARISON

Comparison of Filling Front and Computational Time

The filling in a straight channel is often chosen as the first validation because of its obvious filling pattern. For the comparison of filling front between the new vectorial explicit algorithm and previous algorithm using MINI element, the meshed model and relevant injection parameters are shown in Table 2. In the mesh with MINI tetrahedron elements, there is an extra bubble node in each element, and so it has much more nodes than the mesh of normal tetrahedral elements. However, there is same number of the elements for a model of same discretization. The sliding conditions are imposed on the mold walls. The density of polymer is chosen to be 920 kg [m.sup.-3], viscosity takes a value 250 Pa s, and a velocity equal 1 m [s.sup.-1] is imposed on mold inlet.

TABLE 2. Molded 3D part and injection parameters.

Polymer name           Polypropylene

Polymer temperature    220[degrees]C
Mold temperature       40[degrees]C
In lection velocity    1 m [s.sup.-1]
Viscosity                250 Pa s
Density              920 kg [m.sup.-3]

[TABLE 2 OMITTED]

The comparison of filling patterns obtained by both different algorithms is shown in Fig. 1. The blue bands in postprocessing indicate the positions of filling fronts at different instants when 30%, 60%, and 90% of the mold cavity is filled. It is observed that both algorithms can get similar results in 3D case. This validation is based on the simple models but is effective for evaluating the accuracy of the software.

[FIGURE 1 OMITTED]

To validate the efficiency of the new vectorial algorithm, both the new one and the previous one has been implemented on the MATLAB[C] platform. The geometry model shown in Table 2 is used and meshed from coarse to fine. Using a personal computer with Pentium 4.2 GHz CPU and 256M memory, the models meshed with different number of tetrahedron element have been simulated. The consumed CPU time at each time step for different node number of models is plotted in Fig. 2. It can be observed that the computational cost is about linearly proportional to the degree of freedom number for the new vectorial explicit algorithm, whereas the computational cost is about exponentially proportional to the degree of freedom number for the previous algorithm using MINI elements. Such a feature makes a great difference in computational cost between two algorithms when large industrial problems are considered.

[FIGURE 2 OMITTED]

The major difference between two algorithms locates in the fractional step associated to the incompressibility satisfaction. For the previous explicit algorithm, the global solutions to get pressure fields result in the increase of CPU time about exponentially proportional to the degree of freedom number. The existence of extra bubble node in MINI element results in a computational cost more expensive. In newly developed algorithm, the fields of pressure's gradients are evaluated with the specially designed local and vectorial operations at each time step. There is neither global solution nor the construction of any global matrix, except for the small operations at element level. So the computation can be carried out by merely vectorial operations. This important character makes the CPU time increase with the degree of freedom number much slowly than previous algorithm. Furthermore, the operations in new algorithm are easily to be parallelized in the multicluster systems of high performance computation.

Comparison of the Volume Conservation

The validation of volume conservation during filling process was proposed by Pichelin and Coupez (18). When velocity boundary condition is imposed on the mold inlet, the volume of polymer entering the cavity can be calculated by Eq. 28 at each instant, it can be called theoretical volume of filled polymer.

[Q.sub.t] = [t.summation over [DELTA]t=0] [DELTA]t [[integral] [gamma]' V * [[vector].n dS (28)

where [Q.sub.t] is the volume of polymer filled in the cavity, t corresponds to the injection time, At is the time increment, [bar.V] is the velocity vector imposed on the inlet, and [[vector].n] is the unit normal vector of inlet surface.

At each instant, the field variables corresponding to filling state are obtained by two different algorithms. These results are used to calculate the volume already filled by polymer in the cavity:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (29)

where Vol is the volume filled by the polymer in mold cavity, N are the interpolation functions, and [F.sup.e] are the node values of filling state variable for each elements.

With the simulation model in Table 2, a comparison about the evolution of the filled volumes with respect to time in 3D case is shown in Fig. 3. It shows that the results obtained by both algorithms are globally close to the theoretical obtained ones. The slight differences between the numerical results and the analytical ones mainly resulted from two matters. One matter is the discretization of finite element method. Some small deviation may be issued by the mesh quality of finite element model, which depends on coarseness of the model. Another reason is associated to the initial filling state variable value. To launch the simulation, at least one layer of the elements near the inlet should be set to a value equal 1.0. It causes an initial filled volume in the cavity at the starting stage. However, when large model with fine mesh is simulated, as the cases of industrial application, the difference will decrease.

[FIGURE 3 OMITTED]

Incompressibility Satisfaction

To evaluate the respect of incompressibility in the filled domain, a tapered cavity with increasing section along the filling direction is chosen as an example, shown as Table 3. One can compare the velocity fields obtained from two algorithms and the theoretical result to judge if the incompressibility is maintained in the filled domain in the new algorithm.

TABLE 3. Molded tapered part and injection parameters.

Polymer name           Polypropylene

Polymer temperature        220 C
Mold temperature           40 C
Injection velocity     1 m [s.sup.-1]
Viscosity                250 Pa s
Density              920 kg [m.sup.-3]

[TABLE 3 OMITTED]

The wall condition in mold cavity is supposed to be the sliding one. The density of polymer is chosen to be 4000 kg [m.sup.-3], viscosity takes a value 100 Pa s, and a velocity equal 1 m [s.sup.-1] is imposed on mold inlet. It should be mentioned that these properties and conditions are assigned only for the purpose to evaluate the performance of newly developed software.

Based on simulation results obtained by both algorithms, the curves of axial velocity with respect to axial location are plotted and compared with theoretical result, as shown in Fig. 4. It shows that two methods obtain almost the same results and their results are very close to theoretical one. It proves that incompressibility is well satisfied in the new vectorial algorithm.

[FIGURE 4 OMITTED]

Comparisons with Results of the Previous Algorithm with MINI Element

The injection in a multicavity mold has been chosen to validate the new vectorial algorithm. The mold has three cavities for different components including one wheel part and the corresponding runners are integrated in the same mold. Mesh of the model with triangular elements is shown in Table 4. A pressure equal 16 MPa is imposed on mold inlet, the mold wall is supposed to be the sticking one. A widely used 316L thermal debinding feedstock in Metal Injection Molding is chosen, which is the mixture of 316L stainless steel powder and thermal plastic binder. Density of the feedstock is 5200 kg [m.sup.-3], and the viscous behaviors of the feedstock determined by capillary tests are shown in Fig. 5. The temperature of melted feedstock is 200[degrees]C and the mold temperature is 40[degrees]C.

[FIGURE 5 OMITTED]

TABLE 4. Molded part and injection parameters.

Feedstock name       316L thermal debinding feedstock

Polymer temperature                200 C
Mold temperature                   40 C
Injection pressure                16 MPa
Density                     5200 kg [m.sup.-3]

[TABLE 4 OMITTED]

The comparison of filling patterns at different instants obtained by the previous algorithm with MINI element and the new vectorial algorithm is shown in Fig. 6. It can be observed that the filling patterns obtained by two simulations are generally similar. At the instants corresponding to 12%, 38%, 73%, and 84% of filling ratio, the filling fronts obtained by the new explicit algorithm and the previous one are almost same. The temperature fields obtained by two algorithms at 84% filling ratio is shown in Fig. 7. Two simulation results are almost same except the slight difference of minimum temperature of air in the mold. Despite the sight differences, such a comparison proves that the new developed software performs well for the models, the same as the previous one with MINI element. Moreover, its vectorial feature represents the most important advantages than the performance to simulate the complex models.

[FIGURE 6 OMITTED]

[FIGURE 7 OMITTED]

CONCLUSIONS

A fully explicit vectorial algorithm is proposed and realized for efficient simulation of the mold filling. Under Eulerian description, the incompressibility condition is maintained by a new algorithm without global solution. By virtue of the smoothing procedure, the elements with equal order interpolation for velocity and pressure fields are applied to simplify the industrial applications. The use of MINI elements is no more essential in simulation. Another smoothing procedure is employed to overcome the instability caused by discontinuous mesh filling in the beginning of mold injection. The pressure field can be achieved by a postprocessing for the prescribed instants. As neither global solution nor iteration is required in the new proposed algorithm, the simulation can be performed in a vectorial manner with high efficiency. The computational cost is about proportional to the DOF number of finite element model, instead of the exponential increase with respect to model size in other solution methods. As the numerical operations are performed in a vectorial and uncoupled manner, it can be easily parallelized on a high performance system with multiclusters. Comparison of velocity field obtained by the previous algorithm using MINI element and theoretical solution proves that the feedback regulation strategy in the new algorithm can satisfy well incompressibility. Moreover, the comparison of filling pattern and temperature field results obtained by previous algorithm with MINI element proves the validity of the new algorithm.

NOMENCLATURE

t                       current instant in injection process

X                       spatial vector

[OMEGA]                 mold cavity domain

F                       filling state

[rho]                   density

V                       velocity vector

P                       pressure

[mu]                    viscosity

[epsilon]               shear rate

g gravity               acceleration

[sigma]'                deviatoric Cauchy stress

T                       temperature

[bar.[epsilon]]         equivalent values of shear rate deviators

C                       specific heat

k                       thermal conductivity coefficient

[sigma]': [epsilon]     viscous dissipations

M                       lumped mass matrix

N                       matrix of interpolation function

[F.sup.[sigma]]         effects of viscous diffusion

[F.sup.ext]             external load vector associated to gravity
                        contribution and boundary conditions

[~.G]                   gradient operator

B                       derivative matrix of the interpolation
                        function

[I.sup.e]               constant matrix

[lambda]                factor of regularization

[M.sub.0]               pseudo lumped mass matrix in diagonal form

[K.sup.th]              thermal diffusion term

Q                       heat dissipation term

[q.sup.c]               heat convection term

[Q.sub.t]               volume of polymer filled in the cavity

t                       injection time

[DELTA]t                time increment

[bar.V]                 velocity vector imposed on the inlet

[[vector].n]            unit surface inlet normal vector

Vol                     volume filled by the polymer in mold cavity

REFERENCES

(1.) C.A. Hieber and S.F. Shen, J. Non-Newtonian Fluid Mech., 7, 1 (1980).

(2.) H.H. Chiang, C.A, Hieber, and K.K. Wang, Polym. Eng. Sci., 31(2), 116 (1991).

(3.) G.A.A.V. Haagh and F.N. van de Vosse, Int. .J. Numer. Methods Fluids, 28, 1355 (1998).

(4.) J.-F. Hetu, D.M. Gao, A. Garcia-Rejon, and G. Salloum, Polym. Eng. Sci., 38, 223 (1998).

(5.) F. Ilinca and J.-F. Hetu, Int. Polym. Process., 3. 291 (2001).

(6.) F. Ilinca and J. F. Hetu, Polym. Eng. Sci., 43, 7 (2001).

(7.) R.M. German, Powder Injection Molding-Design and Applications, Pennsylvania University, USA (2003).

(8.) T. Barriere, J.C. Gelin, and B. Liu, .J. Mater. Process. Technol., 143, 636 (2003).

(9.) F. Ilinca, J.F. Hetu, and A. Derdouri, Polym. Eng. Sci., 42, 4 (2002).

(10.) P.M. Gresho, S.T. Chan, R.L. Lee, and CD. Upson, Int.J. Numer. Methods Fluids, 4, 557 (1984).

(11.) R.W. Lewis, H.C. Huang, A.S. Usmani, and J.T. Cross, Int. J. Numer. Methods Eng., 32, 767 (1991).

(12.) T. Barriere, J.C. Gelin, and B. Liu,J. Mater. Process. Technol., 125,518 (2002).

(13.) Z.Q. Cheng, T. Barriere, B.S. Liu, and J.C. Gelin. J. Southwest Jiaotong Univ., 4, 442 (2007).

(14.) T. Hughes, The Finite Element Method, Linear Static and Dynamic Finite Element Analysis. Prentice Hall, Englewood Cliffs (1987).

(15.) R.L. Lee, P.M. Gresho, and R.L. Sani, Int. J. Numer. Methods Eni., 14, 1785 (1979).

(16.) A.S. Usmani, J.T. Cross, and R.W. Lewis. Int.J. Numer. Methods Eng., 35, 787 (1992).

(17.) R.W. Lewis, A.S. Usmani, and T.J. Cross, Int. J. Numer. Methods Fluids. 20, 493 (1995).

(18.) E. Pichelin and T. Coupez, Comput. Methods Appl. Mech. Eng., 163,359 (1998).

Zhiqiang Cheng, (1) Thierry Barriere, (2) Baosheng Liu, (1) Jean-Claude Gelin (2)

(1) Department of Applied Mechanics, Southwest Jiaotong University, 610031 Chengdu, People's Republic of China

(2) FEMTO-ST Institute, ENSMM, Chemin de l'epitaphe, 25030 Besancon, France

Correspondence to: Zhiqiang Cheng: e-mail: [e-mail] Contract grant sponsor: National Natural Science Foundation of China; contract grant number: 10772154: Contract grant sponsor: New Teacher Project of Research Fund for Doctoral Program of High Education of China; contract grant number: 20070613021.

DOI 10.l002/pen.21324

Published online in Wiley InterScience (www.interscience.wiley.com).

[C]2009 Society of Plastics Engineers