Three-dimensional numerical simulation of injection molding filling of optical lens and multiscale geometry using finite element method.INTRODUCTION Injection molding injection molding n. A manufacturing process for forming objects, as of plastic or metal, by heating the molding material to a fluid state and injecting it into a mold. has been an important manufacturing process for the mass production of complex plastic parts and accounts for one-third of all plastics processed. The quality and performance of injection molded parts depend not only on the material, shape, and function of the part design, but also on how the material is processed during molding. By employing scientific-based computer-aided engineering (CAE (1) (Computer-Aided Engineering) Software that analyzes designs which have been created in the computer or that have been created elsewhere and entered into the computer. ) simulation tools, engineers can quickly qualify final design and material selection without physically committing real material and machine time. Despite successful application and wide adoption, the Hele-Shaw flow formulation has its limitations owing to owing to prep. Because of; on account of: I couldn't attend, owing to illness. owing to prep → debido a, por causa de the inherent creeping flow Creeping flow Fluid at very low Reynolds number. In the flow of fluids, a Reynolds number (density · length · velocity/viscosity) describes the relative importance of inertia effects to viscous effects. and thin-wall assumptions. For example, it cannot accurately model three-dimensional flow behaviors within thick and complex components or at the melt fronts (fountain flows), flow junctions, regions where part thickness changes abruptly or where separate melt fronts meet (weld lines), or regions around special part features such as bosses, corners, and/or ribs. Also, a full three-dimensional analysis would be useful in certain applications such as microinjection mi·cro·in·jec·tion n. Injection of minute amounts of a substance into a microscopic structure, such as a single cell. microinjection molded parts, parts with microsurface features or inserts, as well as coinjection and gas-assisted injection molded components because of difficulty and ambiguity in meshing the three-dimensional geometries with thin-shell elements for the 2.5-D simulation [1]. Three-dimensional flow simulation and its applications to a number existing and emerging molding processes have been the active research subjects. Most of the commercial three-dimensional injection molding simulation software Simulation software is based on the process of imitating a real phenomenon with a set of mathematical formulas. It is, essentially, a program that allows the user to observe an operation through simulation without actually running the program. packages [2] are proprietary and being used as "black boxes." Other research-oriented simulation programs tend to be computationally intensive and have been implemented using a specific numerical scheme and finite element See FEA. type. There remains a need to investigate the computational performance and solution accuracy of various numerical schemes and finite element formulations in order to provide useful information to the software companies, the CAE tool users, and the numerical simulation community. In this article, the developments of three-dimensional numerical simulations are presented using different algorithms and finite element types for the filling stage of the injection molding process. The mixed and equal-order formulations are used to solve the Stokes equations, and three different tetrahedral tet·ra·he·dral adj. 1. Of or relating to a tetrahedron. 2. Having four faces. tet elements (namely, Taylor-Hood, MINI, and equal-order) are selected to study their efficiency and accuracy through the non-Newtonian flow tests. A melt tracking algorithm based on the control volume scheme with tetrahedral finite element mesh is used for the three-dimensional simulation of injection molding filling. The equal-order formulation is selected to solve the Stokes equations and the operator splitting method is used to solve the energy equation in the filling simulation of a precision lens. The numerical algorithm is further applied to a molding analysis of a regular sized part with microsurface feature using a two-step, macro-micro filling analysis. MATHEMATICAL MODELING Assuming a creeping flow during the filling in this study, the Stokes equations are the governing equations for the polymer melt flow, which are given as follows: [nabla] u = 0 (1) 0 = - [[nabla].sub.p] + [nabla] (2[eta][dot.[gamma].=]) (2) where [dot.[gamma].=] = [1/2]([nabla]u + ([nabla]u)[.sup.T]). (3) In the above equations, u is the velocity vector, p is the pressure, and [eta] is the viscosity. The heat transfer during the filling is modeled by the energy equation. [rho][c.sub.p]([[partial derivative]T/[partial derivative]t] + u [nabla]T) = [nabla] x (k[nabla]T) + 2[eta][dot.[gamma].=]:[dot.[gamma].=] (4) where [rho] is the density, [c.sub.p] is the specific heat, T is the temperature, and k is the thermal conductivity. The boundary conditions used for solving the Stokes equations in this study are shown below. p = [p.sub.inlet] or u = [u.sub.inlet] on [[GAMMA].sub.inlet] (5) u = 0 on [[GAMMA].sub.wall]. (6) At the melt fronts, the atmospheric pressure atmospheric pressure or barometric pressure Force per unit area exerted by the air above the surface of the Earth. Standard sea-level pressure, by definition, equals 1 atmosphere (atm), or 29.92 in. (760 mm) of mercury, 14.70 lbs per square in., or 101. condition is imposed. The no-slip conditions of Eq. 6 on the mold wall should be imposed only on portions of the cavity filled with polymer melt or unrealistic solutions would occur, as the polymer melt front would never touch the wall. For the energy equation, the initial condition for the entire domain and boundary conditions at the injection gate and mold walls are presented as follows: T = [T.sub.0](x) [for all]x [member of] [OMEGA], t = 0 (7) T = [T.sub.inlet](t) [for all] x [member of] [[GAMMA].sub.inlet], t > 0 (8) T = [T.sub.wall](t) [for all] x [member of] [[GAMMA].sub.wall], t > 0 (9) or q = [h.sub.c]([T.sub.melt] - [T.sub.wall]) [for all] x [member of] [[GAMMA].sub.wall], t > 0 (10) where q is the heat flux at the mold wall and [h.sub.c] is the melt-mold heat transfer coefficient The heat transfer coefficient is used in calculating the convection heat transfer between a moving fluid and a solid in thermodynamics. The heat transfer coefficient is often calculated from the Nusselt number (a dimensionless number). . The boundary condition boundary condition n. Mathematics The set of conditions specified for behavior of the solution to a set of differential equations at the boundary of its domain. at the mold wall is the prescribed temperature or heat flux, and the mold wall boundary condition is imposed only after the wall surface is filled with polymer melt. NUMERICAL IMPLEMENTATION The finite element formulations of the governing equations are described in this section. A control volume method with a tetrahedral element mesh is used as the melt front tracking method. Mixed Formulation In the mixed model, the unknown variables include the velocity and pressure; both types of unknown variables are mixed in a single formulation. Applying the standard Galerkin method In mathematics, in the area of numerical analysis, Galerkin methods are a class of methods for converting an operator problems (such as a differential equation) to a discrete problem. to the Stokes equations and integrating by parts in the momentum equation, the weak formulation Weak formulations have been an important means in the analysis of mathematical equations, since they permit the transfer of concepts of linear algebra to fields like partial differential equations. of the mixed finite element model is constructed as follows [3]: [[integral].sub.[OMEGA]] q[nabla] u d[OMEGA] = 0 (11) [[integral].sub.[OMEGA]] 2[eta][dot.[gamma].=]:[dot.[gamma].=] d[OMEGA] - [[integral].sub.[OMEGA]] p[nabla] v d[OMEGA] = [[integral].sub.[r.sub.c]] ([dot.[tau].=] x n - pn)v ds (12) where [dot.[tau].=] is the deviatoric stress tensor For the stress tensor in classical physics, see the article
[u.sub.i](x) = [N.sub.u.sup.T][~.u.sub.i], p(x) = [N.sub.p.sup.T][~.p] (13) where [~.u.sub.i] (i = 1, 2, 3) and [~.p] are nodal vectors of the velocity components and pressure, respectively. The choice of shape functions for velocity and pressure are different, and is limited by the inf-sup condition (Babuska-Brezzi conditions) that is the mathematical criterion that determines the stability and convergence of a mixed finite element discretization dis·cret·i·za·tion n. The act of making mathematically discrete. [4]. Two types of tetrahedral elements, Taylor-Hood [5] and MINI [6], are chosen to solve the mixed finite element formulations. The quadratic quadratic, mathematical expression of the second degree in one or more unknowns (see polynomial). The general quadratic in one unknown has the form ax2+bx+c, where a, b, and c are constants and x is the variable. and linear tetrahedron tetrahedron: see polyhedron. elements are used for velocity and pressure, respectively, in the Taylor-Hood element. In MINI element, the shape functions for the velocity components are linear interpolation Linear interpolation is a method of curve fitting using linear polynomials. It is heavily employed in mathematics (particularly numerical analysis), and numerous applications including computer graphics. It is a simple form of interpolation. and an introduced higher order bubble function. The shape function for the pressure is linear interpolation. More computational effort is needed for eliminating the internal parameter and the numerical integration In numerical analysis, numerical integration constitutes a broad family of algorithms for calculating the numerical value of a definite integral, and by extension, the term is also sometimes used to describe the numerical solution of differential equations. related to the bubble function compared to the equal-order formulation. Figures 1a-1d show the Taylor-Hood and the MINI elements for velocity and pressure. Note that the mixed velocity-pressure solution does not need to iterate it·er·ate tr.v. it·er·at·ed, it·er·at·ing, it·er·ates To say or perform again; repeat. See Synonyms at repeat. [Latin iter between the velocity and pressure as both are solved simultaneously. For such methods, iterations are needed only because of the nonlinear nature of the equation when non-Newtonian flows are solved, and no iterations would have been needed for a Newtonian fluid. Equal-Order Formulation [FIGURE 1 OMITTED] In the equal-order velocity-pressure formulation [7] for solving the Stokes equations, the same order of discrete form for the pressure and velocity fields is employed. More specifically, the velocity components and pressure are represented using the same linear interpolations of the 4 node tetrahedral element. The velocity and pressure fields are solved by an iterative it·er·a·tive adj. 1. Characterized by or involving repetition, recurrence, reiteration, or repetitiousness. 2. Grammar Frequentative. Noun 1. segregated method. Thus, at any given time, the number of unknowns can be reduced and the memory usage can be restricted. With this implementation, the pressure solution dose not exhibit spurious pressure modes. The standard Galerkin method is used to discretize the momentum equation, and the pressure gradient In atmospheric sciences (meteorology, climatology and related fields), the pressure gradient (typically of air, more generally of any fluid) is a physical quantity that describes in which direction and at what rate the pressure changes the most rapidly around a particular location. term is included on the right hand side in Eq. 14. The discretized momentum equations for each element are assembled to make the global momentum equation shown below (Eq. 15). [[[integral].sub.[OMEGA].sub.c][eta][[partial derivative]N/[[partial derivative][x.sub.j]]], [[[partial derivative][N.sup.T]]/[[partial derivative][x.sub.j]]]dx][~.u.sub.i] + [[[integral].sub.[OMEGA].sub.c] [eta][[[partial derivative]N[partial derivative][N.sup.T]]/[[partial derivative][x.sub.j][partial derivative][x.sub.i]]]dx][~.u.sub.j] = [[integral].sub.[OMEGA].sub.c] N[[partial derivative]p/[[partial derivative][x.sub.i]]]dx (14) [A.sup.u]U = [F.sup.u] or [A.sub.mn.sup.u] [U.sub.n] = [F.sub.m.sup.u]. (15) The pressure equation is derived from the continuity equation by the standard Galerkin method. The pressure equation for each element is assembled to form the following global pressure Eq. 17. More detailed procedures, including boundary conditions, can be found in Ref. 7. After solving the pressure equation, the solutions are used to update the velocity. [[[integral].sub.[OMEGA].sub.c] [[partial derivative]N/[[partial derivative][x.sub.j]]]([N.sup.T][~.K.sub.i])[[[partial derivative][N.sup.T]]/[[partial derivative][x.sub.j]]]dx][~.p] = [[[integral].sub.[OMEGA].sub.c][[partial derivative]N/[[partial derivative][x.sub.i]]][N.sup.T]dx][^.u] + [[line integral around closed path].sub.[GAMMA].sub.c] N[u.sub.i]ds (16) [A.sup.p]P = [F.sup.P] (17) where [^.U.sub.m] = - [summation summation n. the final argument of an attorney at the close of a trial in which he/she attempts to convince the judge and/or jury of the virtues of the client's case. (See: closing argument) over (n[not equal to]m)][[[A.sub.mn.sup.u][U.sub.n]]/[A.sub.mn.sup.u]] (18) [K.sub.m] = [^.K.sub.m]/[A.sub.mm.sup.u]([^.K] = [summation over c][[integral].sub.[OMEGA].sub.c] N dx). (19) Operator Splitting Method The low thermal diffusivity In heat transfer analysis, thermal diffusivity (symbol:  ) is the ratio of thermal conductivity to volumetric heat capacity.n. 1. The transfer of a property of the atmosphere, such as heat, cold, or humidity, by the horizontal movement of an air mass: equation and a diffusion equation The diffusion equation is a partial differential equation which describes density fluctuations in a material undergoing diffusion. It is also used to describe processes exhibiting diffusive-like behaviour, for instance the 'diffusion' of alleles in a population in population as follows [8] [[T* - [T.sub.0]]/[DELTA]t] + u[nabla]T* = 0 (20) [rho][c.sub.p]([[T - T*]/[DELTA]T]) = [nabla] (k[nabla]T) + 2[eta][dot.[gamma].=]:[dot.[gamma].=] (21) where [T.sub.0] is temperature at the previous time step, T* is an intermediate solution of the advection equation, T is the temperature at the present time step, and [DELTA]t is the time step. This method solves the advection equation first and then solves the diffusion equation. A flux-corrected transport (FCT FCT Faculdade de Ciências e Tecnologia (Portuguese University) FCT Fundamentals of Computation Theory FCT Fundação para a Ciência e a Tecnologia (Portuguese Science and Technology Foundation) ) approach [9] is used to solve the advection Eq. 20 with an implicit time discretization. The solution is determined by the proper combination of both solutions of a higher-order scheme and a lower-order scheme. A higher-order scheme of the first solution provides accuracy while a lower-order scheme of the second solution contains artificial viscosity to reduce nonphysical oscillations [8]. The details of this limiting procedure for the combination of the both solutions are described in the Ref. 9. The Galerkin/gradient-least-squares (GGLS) formulation [10] is used to solve the diffusion Eq. 21. This formulation is an effective method for problems with a material that has a low conductivity and thin thermal boundary layers. Control Volume Method with Tetrahedral Element Mesh This approach is one of the volume methods that use a fixed reference for tracking free moving boundaries at the polymer melt fronts. It is a new approach that extends the conventional two-dimensional flow analysis network (FAN) method to three-dimensional filling simulations. In this method, a single control volume comprised several subelement portions of the finite elements associated with a particular node. Each element is divided into four subelements as shown in Figure le and each subelement is a portion of the control volume of a node that each subelement shares. The control volume of the node i is defined as [V.sub.cv.sup.i] = [ni.summation over j][([V.sub.sub_element.sup.i].sub.j] (22) where [V.sub.sub_element.sup.i] is the volume of a subelement that has node i and [n.sub.i] is the total number of subelements associated with node i. For the control volume of a particular node, the value of the fill factor is given and the fill factor for node i in the mesh is defined as follows: [f.sub.i] = [V.sub.ff.sup.i]/[V.sub.cv.sup.i] (23) where [V.sub.cv] is the volume of the control volume and [V.sub.ff] is the volume of filled fluid in the control volume. In this method, the fin factor, f, is defined as f = 1 when the control volume is filled with polymer and f = 0 when the control volume is empty. The fill factor, f, has values between 1 and 0 in the partially filled control volumes of the melt fronts in the cavity. To update the fill factor of all partially filled control volumes, the flow rate of the polymer melt into the control volume should be calculated. The calculation of flow rate is done for each subelement. The summation of flow rates for all of the subelements is the total flow rate of a control volume. When calculating the flow rate for a subelement, only the flow rates from the triangular surfaces that contact fully filled control volumes are considered. The flow rate through a contacting triangular surface is evaluated by a surface integral of normal velocity over the triangular area. On the basis of the flow rates of all partially filled control volumes, the shortest time to fully fill one of the partially filled control volumes is selected as the time increment to update the fill factors. After one of the partially filled control volumes is fully filled by updating the fill factors, the computational domain is also updated. All elements that have the node corresponding to the fully filled control volume are included into the computational domain and the newly included nodes become the new partially filled nodes. The melt fronts can be determined as iso-value surfaces of fill factors defined at nodes in the finite element mesh. In this study, the iso-value surfaces of f = 0.5 are considered as melt fronts. The iso-value surfaces are obtained using the nodal values of fill factors and linear interpolation functions for each finite element in the mesh. Generally, three-dimensional mold filling simulations require solving the polymer melt flow equations and energy equation while tracking the melt fronts at every time step until the cavity is filled completely. In the presentation implementation, the number [N.sub.f] is the number of fully filled nodes, which is counted for each iteration, and [N.sub.c] is the number of nodes that need to be filled before solving the Stokes equations, which is set as an input data. The minimum for [N.sub.c] is the number of nodes that are fully filled when at least one element is added to the computational domain or the boundary conditions are changed. To speed up the calculation of the filling simulation, more nodes can be filled than the minimum [N.sub.c]. The constant time step for the energy equation is given from an input file and, using this time step, solving the energy equation is iterated with the melt tracking algorithm until the number of filled nodes reaches [N.sub.c]. The integration time step for the energy equation is set as a constant based on the Refs. 8 and 9 while the time step for filling varies based on the melt fronts tracking algorithm. Detailed description on the algorithm, the flow chart, and the determination of the time step for filling simulation used in this study can be found in [11]. Memory Management The global stiffness matrix resulting from the finite element method is typically a large sparse, banded symmetric, or nonsymmetric matrix system, especially in three-dimensional cases. In this study, the conjugate conjugate /con·ju·gate/ (kon´jdbobr-gat) 1. paired, or equally coupled; working in unison. 2. a conjugate diameter of the pelvic inlet; used alone usually to denote the true conjugate diameter; see gradient (CG) solver is used to solve the stiffness matrices. In solving the large sparse system, extensive CPU time The amount of time it takes for the CPU to execute a set of instructions and generally excludes the waiting time for input and output. CPU time - processor time and memory storage requirements become important issues. A couple of large sparse matrix In the mathematical subfield of numerical analysis a sparse matrix is a matrix populated primarily with zeros. Conceptually, sparsity corresponds to systems which are loosely coupled. solvers and storage schemes for finite element method analysis have been developed as reported in [3, 12]. In this study, a modified procedure of the conventional assembling method for the global stiffness matrix was proposed to employ the CG solver without using a huge amount of storage space. This memory management procedure selects a small block of coefficient within the banded matrix and moves them down as shown in Fig. 2. When the block is created, the global matrix information within the block region is extracted, including the boundary conditions. Before moving down the block, all non-zero values are stored in another array with their row and column information. Based on this simple yet effective implementation, the CPU time can be reduced considerably without using the computer virtual memory. For example, using a computer with available physical memory of 270 MB, the assembling time for a 630 MB size banded symmetric matrix symmetric matrix n. Mathematics A matrix that is its own transpose. was about 40 sec with 200 MB block size while it took 17 min with the virtual memory. [FIGURE 2 OMITTED] NUMERICAL VERIFICATIONS AND VALIDATION RESULTS Non-Newtonian Flow Tests Non-Newtonian flow tests in simple geometries with known analytical solutions were conducted as the verification and comparison studies for the present numerical implementation with three different tetrahedral elements (i.e., Taylor-Hood, MINI, and equal-order). The geometries with the fluid flowing within and the boundary conditions for the flow tests are shown in Fig. 3. The power-law model was used as the viscosity model in non-Newtonian, isothermal flow Isothermal flow is a model of a fluid flow, which remains in the same temperature. In this model the temperature is constant while the stagnation temperature is changing. The change in stagation temperature occur because the temperature is constant but the velocity increasing. tests as the analytical solution is readily available for comparison and verification. [eta] = m[dot.[gamma].sup.n-1] (24) where [dot.[gamma]] is strain rate. The material constants used in these tests were based on those of a high-density polyethylene high-density polyethylene n. Abbr. HDPE A strong, relatively opaque form of polyethylene having a dense structure with few side branches off the main carbon backbone. (HDPE HDPE abbr. high-density polyethylene ) and the values of m and n in the model are 2.374 X [10.sup.4] (Pa [s.sup.n]) and 0.424, respectively. At the inlet and exit surfaces, pressure boundary conditions of [p.sub.i] (= [10.sup.5] Pa) and [p.sub.o] (= 0 Pa) were applied to both the plate and cylindrical geometries. The numbers of elements and nodes in the finite element mesh for the plate and cylindrical models are shown in Tables 1. For the numerical tests, a computer with Intel ZEON (3.0 GHz) processor and 1.5 GB memory was used and the usage of physical memory was allowed up to 1.0 GB to allocate the variables of the global stiffness system matrix in the tests. The relaxation factor was set at 0.5, and the calculation was repeated with 40 iterations (regardless whether the convergence criterion is met or not) in all tests of the three element sets to compare their convergence characteristics. The error criterion for the conjugate gradient solver is set to be 1.0 X [10.sup.-6]. [FIGURE 3 OMITTED] The velocity and pressure distributions and velocity profiles at the locations of inlet (z = 0 mm), middle (z = 30 mm), and exit (z = 60 mm) using the Taylor-Hood element are shown in Figs. 4 and 5 for the plate and cylindrical models with the analytical solutions. The predicted pressure profiles showed a very uniform pressure gradient along the flow direction. The predicted velocity profiles at all three different locations were almost identical and quite accurate comparing the analytical solutions for both plate and cylindrical models. The errors of the velocity profiles after 40 iterations using the MINI element and the equal-order element are compared with the errors of the Taylor-Hood element in Table 2. The errors of the MINI element and the equal-order element were about 4% near y = 0 mm and 3% range near the center of r = 0 mm in plate and cylindrical models, respectively, and the errors increased near the wall regions at y = [+ or -]10 mm and r = 10 mm. These trends in error were similar to that of the Taylor-Hood element case, but the results of Taylor-Hood element were more accurate than those of MINI and equal-order elements, which had 1% range errors. In the wall regions at y = [+ or -]10 mm (plate) or r = 10 mm (cylinder), where the velocity variations were much larger and the velocities approached to zero rapidly so that the errors became large. It is expected that fine meshes with increased number of elements or an adaptive meshes with better resolution near the mold wall can help reduce the errors in these mold wall regions. This table shows the CPU CPU in full central processing unit Principal component of a digital computer, composed of a control unit, an instruction-decoding unit, and an arithmetic-logic unit. times to convergence and the average errors of the velocity at all nodal points after 40 iterations of each element set. The CPU times to convergence were defined as times when the errors were within the 1% range of the errors in 40 iterations. Although the most accurate results were obtained after 40 iterations, Taylor-Hood element was too expensive in CPU times to be used for solving the governing equations in three-dimensional filling simulations. The equal-order element used the smallest times to convergence in 40 iterations. Although the node and element numbers for the MINI and equal-order elements were the same, the CPU times of the MINI element were longer than those of the equal-order element. The average errors of the MINI and equal-order elements were almost at the same level. This means that the equal-order formulation is more efficient method than the mixed formulation with the MINI element in solving the Stokes equations. The shorter CPU time with the equal-order element was because the equal-order formulation used a separate approach for the velocity and pressure equations in solving the Stokes equations, and the matrix sizes in the formulation and CPU times used for assembling the matrices were smaller than those of the mixed formulation with the MINI element. Based on these numerical tests and results comparisons, it was found that the Taylor-Hood element gave the most accurate results but with the longest CPU time. The results from the MINI element and the equal-order element are comparable but the latter required less CPU time. The equal-order formulation is three to nine times faster than MINI element in convergence rate as seen in Table 2. Detailed results of numerical solution convergence and the rate of convergence In numerical analysis (a branch of mathematics), the speed at which a convergent sequence approaches its limit is called the rate of convergence. Although strictly speaking, a limit does not give information about any finite first part of the sequence, this concept is of practical can be found in [11]. [FIGURE 4 OMITTED] [FIGURE 5 OMITTED] Three-Dimensional Filling Simulation for an Optical Lens Part The three-dimensional filling was performed for an industrial precision optical lens part. The equal-order formulation was selected to solve the Stokes equations, and the energy equation was solved using the operator splitting method in the filling simulation. Experiments. Experiments for the injection molded lens part were set up to produce several short shot or full shot lens parts to validate the predicted results. The experiments for the injection molded optical lens parts were performed at 3M Precision Optics. The parts were molded in a four-cavity mold, and the parts were located at the end of four runners. The material used in the experiments was PMMA PMMA polymethyl methacrylate. Plexiglas V-825 from Rohm & Hass. The outside diameter Outside diameter is the diameter of the addendum (tip) circle. In a bevel gear it is the diameter of the crown circle. In a throated wormgear it is the maximum diameter of the blank. The term applies to external gears.1 Notes 1. of the lens part is 96.19 mm, and the height of part at the center is 19.87 mm. The largest thickness of the lens part is at the outer rim of the lens and is 10.50 mm while the thickness at the center is 6.00 mm. The weight of each lens part is 69.8 g. The melt and mold wall temperature used in the injection molding experiments were 505.37 and 349.3 K, respectively. In these experiments, short shots of one fourth, one half, and three fourth of the fully filled volume were collected. As a common practice used in industry, the shapes of short shot parts were used to portray the filling patterns during the filling and compared with the simulation results for the validations. Material Properties. The rheological rhe·ol·o·gy n. The study of the deformation and flow of matter. rhe  o·log behavior of the PMMA was
described by a 7-constant modified-cross model with WLF WLF Washington Legal FoundationWLF Wallis and Futuna (ISO Country code) WLF Waist Level Finder (camera viewfinder type) WLF Viva La Figa (MotoGP motorcycle races) equation, which is represented as follows: [eta]([dot.[gamma]],T,p) = [[[eta].sub.0](T,p)]/[1 + ([[[eta].sub.0][dot.[gamma]]]/[tau]*)[.sup.1-n]] (25) [[eta].sub.0](T,p) = [D.sub.1] exp exp abbr. 1. exponent 2. exponential (-[[[A.sub.1](T - T*)]/[[A.sub.2] + (T - T*)]] (26) T* (p) = [D.sub.2] + [D.sub.3]p (27) [A.sub.2] = [~.A.sub.2] + [D.sub.3]p. (28) For PMMA, the viscosity parameter values of n, [tau]*, [A.sub.1], [~.A.sub.2], [D.sub.1], [D.sub.2], and [D.sub.3] are 0.26359, 9.9338 X [10.sup.4] Pa, 42.565, 51.6 K, 5.86 X [10.sup.16] Pa s, 377.15 K, and 0 K [Pa.sup.-1]. The material properties of [rho], k, and [c.sub.p] are 1185.0 kg/[m.sup.3], 0.21 W/m K, and 2300 J/kg K. [FIGURE 6 OMITTED] Model and Simulation. The finite element mesh and boundary conditions for a half of the optical lens (due to symmetry) are shown in Fig. 6. The finite element mesh was generated using 4-node tetrahedral elements. The mesh had 90,352 tetrahedral elements and 17,355 nodes. There were 4-5 element layers across the lens thickness on the symmetric plane. Based on the process conditions of the experiments, a uniform inlet velocity of 120 mm/s was specified so that the fill time for lens part would be 4.9 sec. The inlet temperature of 505.37 K was imposed at the gate, which was the melt temperature used in the experiments. The heat transfer coefficient of 4,000 W/([m.sup.2] K) was applied at the mold wall surfaces. The velocity vectors and pressure distributions near the end of the filling (at t = 4.8 sec) are shown in Fig. 7. Each surface represents a uniform pressure plane of 0, 1, 2, 3, and 4 MPa, respectively. It is seen that the pressure gradient is higher in the center area of the part than that of the thicker edge area. The figure also shows that the material moves faster in the edge regions of the lens part, which are thicker than the center regions. This race-tracking phenomenon in the edge region resulted in a special reverse "Y" weld line pattern near the last-to-fill region. The predicted temperature distribution (not shown) on the symmetric plane at various fill times revealed that the temperature at the interface of the melt and mold decreased below the glass-transition temperature but the thermal boundary layer boundary layer In fluid mechanics, a thin layer of flowing gas or liquid in contact with a surface (e.g., of an airplane wing or the inside of a pipe). The fluid in the boundary layer is subjected to shear forces. was so thin that the boundary layer was within one element layer because of the low thermal conductivity of the polymer melt. Within the part, local high temperature areas due to shearing heating was found. [FIGURE 7 OMITTED] [FIGURE 8 OMITTED] Experimental Validations. Figure 8 shows the short shot parts from the experiments and predicted filling patterns of the corresponding positions from the simulation. Comparing with the short shot parts, the filling patterns of simulation results for one fourth, one half, and three fourth filled parts showed very good agreements. Slight differences of melt advancements were found between experimental and simulation results. The discrepancy is due to volumetric volumetric /vol·u·met·ric/ (vol?u-met´rik) pertaining to or accompanied by measurement in volumes. vol·u·met·ric adj. Of or relating to measurement by volume. shrinkage of the actual molded parts, numerical truncation errors, and modeling of the material properties, process conditions, and geometry and mesh. An unusual pattern of fainted weld lines resembling a reversed "Y" was found in a full shot lens part and was highlighted with black lines in Fig. 9. The weld lines were predicted in the numerical results, and Fig. 9 shows the formation processes at different fill times before the end of filling. [FIGURE 9 OMITTED] It has been suggested that a meeting angle at the advancing melt front ([theta Theta A measure of the rate of decline in the value of an option due to the passage of time. Theta can also be referred to as the time decay on the value of an option. If everything is held constant, then the option will lose value as time moves closer to the maturity of the option. ]) (cf. Fig. 9) of less than 135[degrees] is thought to have a high possibility of forming a weld line on the part surface [13]. Based on simulation prediction, a slanted weld line segment first formed as the melt from the thinner center region and the met fronts from the thicker edge regions met. Next, two vertical fronts from the left and right sides of the lens met and formed a vertical segment of the weld line. Using the simulation, a number of process condition modifications were evaluated to reduce the degree of race-tracking and eliminate/reduce the weld line by increasing the meeting angle. Figure 10 shows more desirable filling patterns by using the higher mold wall temperature of 399.3 K. Although a higher mold wall temperature will potentially increase the cycle time, it is not uncommon in precision lens molding to temporarily and locally heat up mold surfaces prior to filling to eliminate residual stresses and birefringence Birefringence The splitting which a wavefront experiences when a wave disturbance is propagated in an anisotropic material; also called double refraction. In anisotropic substances the velocity of a wave is a function of displacement direction. without causing a significant increase in cooling time (Law) such a lapse of time as ought, taking all the circumstances of the case in view, to produce a subsiding of passion previously provoked. - Wharton. See also: Cooling . FILLING SIMULATION FOR A PART WITH A MICROSURFACE FEATURE The plastic microcomponents have great potential in fabricating the miniaturized devices of electro/optical, mechanical, and biomedical bi·o·med·i·cal adj. 1. Of or relating to biomedicine. 2. Of, relating to, or involving biological, medical, and physical sciences. applications. In addition, many macro-sized devices in those applications also contain high aspect ratio (depth-to-width) microsurface features. The accurate replication of the microsurface features is important in the fabrication fabrication (fab´rikā´sh  n),n the construction or making of a restoration. process. The material properties, process conditions, and geometry of the microsurface features are key factors that influence the molding quality of the micro-surface feature [14]. In this section, the filling of a plate part with a microsurface feature is analyzed through the three-dimensional simulations. [FIGURE 10 OMITTED] Two-Step Macro-Micro Filling Analysis The plastic parts in microinjection molding are either micrometer-sized parts or regular-sized parts with micro features, with the latter involving multiscales in geometries and mesh sizes. In three-dimensional case, very fine elements are needed to mesh the micro features, which would result in an excessive number of elements if the macro-scaled substrate were to be discretized with the same or adaptive mesh size. A 2.5-D/3-D hybrid numerical technique has been used to reduce the elements in the mesh and computation time In computational complexity theory, computation time is a measure of how many steps are used by some abstract machine in a particular computation. For any given model of abstract machine, the computation time used by that abstract machine is a computational resource which can be [15]. In this technique, three-dimensional simulation was used around the micro feature whereas a 2.5-D simulation was used for the base plate. In a separate, parallel study, a commercial injection molding software, C-MOLD, was used to analyze the filling and postfilling of a high aspect ratio micro feature on a thick substrate. 2.5-D simulations were used for both the micro feature and the substrate but a two-step approach was proposed to avoid direct modeling of both the substrate and the micro feature in the simulation [14]. In this study, a three-dimensional, two-step macro-micro filling analysis is introduced to simulate the filling of the part with a microsurface feature. First, a filling analysis is performed using the macro-scale geometry of the part. Since the micro feature hardly affects the filling of the macro-level part, it can be ignored in the macro-level analysis. Next, a micro-level filling analysis on the micro-size geometry is performed using the results from the macro-level analysis as the input and boundary conditions. The main goal of this two-step approach is to avoid an excessive number of finite elements in the mesh dictated by the micro feature. Also, more accurate and detailed information can be obtained from the simulation of the separated micro-level analysis with smaller time step than that of macro-level analysis. As the size of the micro feature decreases, the temperature variation at the inlet of the micro feature can be assumed to be very small. Thus, the effect of possible three-dimensional cooling around the micro feature inlet on the melt temperature distribution is ignored in this study. Figure 11 shows the geometry of the micro-sized surface feature on the base plate based on Ref. 15 and the procedure of a two-step filling analysis. The base plate is a 2-mm thick and 41-mm long rectangular plate, and the material is injected into the cavity through a gate with constant flow rates. The micro feature is a rectangular plate with an aspect ratio of 4 (thickness: 200 [micro]m, length 800 [micro]m). The capillary number In fluid dynamics, the capillary number represents the relative effect of viscous forces versus surface tension acting across an interface between a liquid and a gas, or between two immiscible liquids. It is defined as vis·cous adj. 1. Having relatively high resistance to flow. 2. Viscid. stress to surface tension in a free surface flow, is defined as [FIGURE 11 OMITTED] [C.sub.a] = [eta]U/[sigma] (29) where [eta], [sigma], and U are viscosity, surface tension, and characteristic velocity, respectively. On the basis of the numerical results with inlet velocity, [U.sub.o] = 25 mm/sec and [h.sub.c] = 1000 W/[m.sup.2] K for the base plate, the average local velocity, pressure, temperature, and shear rate Shear rate is a measure of the rate of shear deformation:  For the simple shear case, it is just a gradient of velocity in a flowing material. were assumed to be 1.577 mm/s, 0.1 MPa, 480 K, and 50 [sec.sup.-1], respectively, during the filling in the micro feature. The capillary number was estimated to be over 100, and hence the surface tension was negligible when compared with viscous stress. But, as the temperature and shear rate increase, the capillary number tends to decrease. In the micro feature, there is high possibility of the wall-slip because of the fast build-up of the shear stress shear stress n. See shear. shear stress A form of stress that subjects an object to which force is applied to skew, tending to cause shear strain. at the mold wall surfaces. If the wall-slip condition were imposed on the mold surfaces, the filled heights would be higher with the increased average Mow velocity in the micro feature [15]. At present, effects such as surface tension, surface roughness, viscoelastic Adj. 1. viscoelastic - having viscous as well as elastic properties natural philosophy, physics - the science of matter and energy and their interactions; "his favorite subject was physics" effect, cooling by heat loss at the melt front, etc. for the micro filling were ignored. Nevertheless, it may be necessary to include these effects in future studies for a more accurate prediction for the flow behavior within the microsurface feature. Numerical Results for the Base Plate Since the flow in the base plate is essentially one-dimensional, a small strip of the plate along the length, which has 2-mm width, 2-mm thickness, and 41-mm length, was used as the computation domain instead of the whole plate cavity to reduce the computational time. The uniform flow velocities ([U.sub.o]) in the base plate were 3.125, 6.25, 12.5, 25, 50, and 75 mm/sec. Inlet melt temperature ([T.sub.o]) and mold temperature ([T.sub.w]) in base plate were constant as 500.15 K and 298.15 K, respectively. Constant heat transfer coefficients with values of 500, 1000, 2000, and 3000 W/[m.sup.2] K were used. The numbers of elements and nodes were 26,059 and 5,604 in the mesh. The mesh had about 6 element layers in thickness direction. The material used in the numerical simulation was PMMA, PL-150 from Plaskolite. The glass-transition temperature ([T.sub.g]) is 373.15 K, and the material properties of [rho], k, and [c.sub.p] are 1060.0 kg/[m.sup.3], 0.15 W/m K, and 2050.0 J/kg K. The viscosity parameter values of n, [tau]*, [A.sub.1], [~.A.sub.2], [D.sub.1], [D.sub.2], and [D.sub.3] in the Cross-WLF model are 0.3791, 2.01 X [10.sup.4] Pa, 51.44, 51.6 K, 1.37 X [10.sup.20] Pa s, 337.1 K, and 0 K [Pa.sup.-1], respectively [15]. The predicted temperature and pressure histories at the location of micro feature with different process conditions were recorded, which were then used as input data for the filling simulation of the micro feature. The inlet pressure and temperature profiles for the micro feature versus filled volume in the base plate with different velocities are shown in Fig. 12. The pressure increased from zero at near 50% filled volume until the base plate was fully filled. The pressure increase was larger as the flow speed in the base plate decreased, reflecting the effect of cooling. Meanwhile, the temperature near the micro feature decreased gradually during the filling of the base plate. This figure shows that the material near micro feature was colder with decreasing velocity. The inlet pressures for micro feature versus fill time at the end of filling are shown in Fig. 13a. The longer fill time and higher heat transfer coefficient (more effective cooling at the melt-mold interface) resulted in higher pressure at the end of the filling. The inlet temperature at the end of filling for micro feature versus fill time is shown in Fig. 13b. At the end of the filling, the inlet temperatures at the micro feature were cooled for below [T.sub.o] and in some cases even below [T.sub.g]. This figure clearly shows that polymer melt entering the micro feature was colder with longer fill times and higher heat transfer coefficients. To continuously fill the micro feature, the temperature of the material should be higher than [T.sub.g]. If the material inside the micro feature is fully solidified, the material cannot fill the micro feature anymore even though the inlet temperature is higher than [T.sub.g]. One can see that the flow velocity In fluid dynamics the flow velocity, or velocity field, of a fluid is a vector field which is used to mathematically describe the motion of the fluid. Definition The flow velocity of a fluid is a vector field [FIGURE 12 OMITTED] Numerical Results for the Micro-Surface Feature For numerical analysis numerical analysis Branch of applied mathematics that studies methods for solving complicated equations using arithmetic operations, often so complex that they require a computer, to approximate the processes of analysis (i.e., calculus). of the filling for the micro feature, 28551 tetrahedral elements and 5636 nodes were used, and the mesh had about 6 element layers in the thickness direction. The variable pressure and temperature, which were obtained from macro-level filling simulation of the base plate, were imposed on the micro feature inlet surface. The ambient mold wall temperature was 298.15 K, and the same constant heat transfer coefficients as the base plate were used. With [U.sub.o] = 75 mm/s and [h.sub.c] = 1000 W/[m.sup.2] K for the base plate, the melt front advancements in the micro feature during the filling are shown in Fig. 14. Initially, the height of the melt front position increased rapidly and then the advancement slowly came to a stop in the middle of the micro feature. The 800 [micro]m height micro feature could not completely filled after the filling was finished in 0.5467 sec. Figure 15 shows the temperature distributions in the cross-section at the center of the micro cavity during filling. At the end of the filling, the temperature of the material was over [T.sub.g] and still high enough to flow inside the micro feature. However, the inlet pressure was not large enough to push the material into the cavity any further and eventually the material inside the micro feature solidified completely. [FIGURE 13 OMITTED] [FIGURE 14 OMITTED] For the validation of the numerical predictions, the filled lengths in micro feature versus the flow velocity with different heat transfer coefficients were compared with the experimental data in Fig. 16. Experimental data were obtained from the reference [15]. Based on the comparison, the three-dimensional, two-step filling analysis predicted well the filling in the micro feature and the predicted filled lengths in the micro feature correlated well with the experimental data. When the flow velocity was less than 20 mm/sec, the experimental data was coincident well with the numerical results of [h.sub.c] = 500 W/[m.sup.2] K. However, if the flow velocity is higher than 20 mm/sec, the experimental data closely followed the numerical results of [h.sub.c] = 1000 W/[m.sup.2] K. This suggests that the heat transfer coefficient may depend on the flow velocity, which influences the pressure and temperature inside the micro feature. It is clearly seen that the flow velocity in the base plate and heat transfer coefficient affected the filled lengths in micro feature. A higher injection speed reduced the fill time of the base plate. It also reduced the frozen layer and prevented premature freezing at the entrance of the micro feature. In addition, the temperature of the material entering the micro feature is also higher so that the micro feature could be filled more. A small heat coefficient resulted in higher filled length in the micro feature and the change of the heat transfer coefficient resulted in a large difference in the filled lengths. Since the material in the micro feature and the base plate was warmer with a smaller coefficient value, the temperature of the material in the micro feature was higher. [FIGURE 15 OMITTED] [FIGURE 16 OMITTED] CONCLUSIONS The three-dimensional numerical simulations for the filling stage of the injection molding process were studied using different algorithms and finite element types. The equal-order model was more efficient than the mixed model of the MINI element in solving the Stokes equations because of the separate approach in solving the velocity and pressure equations. A simple memory management procedure was introduced to deal with the large sparse matrix system and reduced the CPU time without using a huge amount of storage space. Through the experimental validation using a three-dimensional industrial precision optical lens, it could be seen that the numerical implementation worked well in predicting the melt front advancements using the control volume scheme with tetrahedral finite element mesh. The numerical algorithm for filling was applied to the filling analysis of a rectangular part with a micro feature with a three-dimensional, two-step macro-micro filling approach. The predicted filled lengths in the micro feature correlated well with the experimental data and suggested that variable heat transfer coefficients should be used for a more accurate prediction for parts with micro-surface features. ACKNOWLEDGMENTS The authors thank Roy Auerbach, James P. Ruth, and Kathleen Allen of 3M for their useful discussion, assistance, and the experimental data. REFERENCES 1. S.W. Kim and L.S. Turng, Modell. Simul simul /sim·ul/ (sim´ul) [L.] at the same time as. . Mater. Sci. Eng., 12, S151 (2004). 2. Moldflow Corporation, http://www.moldflow.com/Products/MPI/modules/3d.asp and Coretech System, http://www.moldex3d.com/ 3. J.N. Reddy and D.K. Gartling, The Finite Element Method in Heat Transfer and Fluid Dynamics fluid dynamics n. (used with a sing. verb) The branch of applied science that is concerned with the movement of gases and liquids. , CRC (Cyclical Redundancy Checking) An error checking technique used to ensure the accuracy of transmitting digital data. The transmitted messages are divided into predetermined lengths which, used as dividends, are divided by a fixed divisor. Press, Boca Raton Boca Raton (bō`kə rətōn`), city (1990 pop. 61,492), Palm Beach co., SE Fla., on the Atlantic; inc. 1925. Boca Raton is a popular resort and retirement community that experienced significant industrial development in the 1970s and 80s. , FL (1987). 4. K.J. Bathe, Finite Element Procedures, Prentice-Hall, Englewood Cliffs, NJ (1996). 5. P. Hood and C. Taylor, "Navier-Stokes Equations Using Mixed Interpolation interpolation In mathematics, estimation of a value between two known data points. A simple example is calculating the mean (see mean, median, and mode) of two population counts made 10 years apart to estimate the population in the fifth year. ," in Finite Element Methods in Flow Problems. J.T. Oden, R.H. Zienkiewicz, R.H. Gallagher, and C. Taylor, editors, UAH press, University of Alabama, Huntsville, AL, 121 (1974). 6. D.N. Arnold, F. Brezzi, and M. Fortin, Calcolo, 21, 337 (1984). 7. J.G. Rice and R.J. Sehnipke, Comp. Methods Appl. Mech. Eng., 58, 135 (1986). 8. F. Ilinca and J.F. Hetu, Int. J. Numer. Meth. Eng., 53, 2003 (2002). 9. R. Lohner, K. Morgan, J. Peraire, and M. Vahdati, Int. J. Numer. Meth. Fluid., 7, 1093 (1987). 10. L.P. Franca, E.G. Dutra do Carmo, Meth. Appl. Mech. Eng., 74, 41 (1989). 11. S.W. Kim, "Three-Dimensional Simulations for the Filling Stage of the Polymer Injection Molding Process using Finite Element Method", Ph.D. Thesis, University of Wisconsin-Madison “University of Wisconsin” redirects here. For other uses, see University of Wisconsin (disambiguation). A public, land-grant institution, UW-Madison offers a wide spectrum of liberal arts studies, professional programs, and student activities. , Madison, WI (2005). 12. D.T. Nguyen, Parallel-Vector Equation Solvers for Finite Element Engineering Applications. Kluwer, Boston (2002). 13. L.S. Turng, B. Nagel, and N.S. Liu, C-MOLD Design Guide, C-MOLD. Ithaca, NY (1994-2000). 14. D. Yao and B. Kim, J. Injection Molding Technol., 6, 11 (2002). 15. L. Yu, L.J. Lee, and K.W. Koelling, Polym. Eng. Sci., 44, 1866 (2004). Sang-Woo Kim, Lih-Sheng Turng Polymer Engineering Center, Department of Mechanical Engineering, University of Wisconsin-Madison, Madison, Wisconsin Madison is the capital of the U.S. state of Wisconsin and the county seat of Dane County. It is also home to the University of Wisconsin–Madison. The 2006 population estimate of Madison was 223,389, making it the second largest city in Wisconsin, after Milwaukee, and 53706 Correspondence to: L.-S. Turng: e-mail: turng@engr.wisc.edu Contract grant sponsor: 3M Precision Optics, Inc.; contract grant sponsor: National Science Foundation; contract grant number: DMI-0323509.
TABLE 1. The numbers of elements and nodes for three tetrahedral
elements in the plate and cylindrical models.
Node Node
Element (velocity) (pressure)
Plate model
Taylor-Hood 10,043 14,804 2031
MINI 10,043 2031 2031
Equal-order 10,043 2031 2031
Cylinder model
Taylor-Hood 10,003 14,533 1972
MINI 10,003 1972 1972
Equal-order 10,003 1972 1972
TABLE 2. The CPU times to convergence and average errors of velocity at
all nodes in 40 iterations for three tetrahedral elements.
Plate Cylinder
Avg.
CPU time error CPU time Avg. error
Element (sec) (%) (sec) (%)
Taylor-Hood 1.61 x [10.sup.4] 1.61 6.23 x [10.sup.3] 2.20
MINI 375.22 5.67 180.27 5.03
Equal-order 39.13 5.72 55.98 4.91
|
|
||||||||||||||||
Printer friendly
Cite/link
Email
Feedback
Reader Opinion