Thermal flow simulation for casting/molding processes

[0001] The present application relates to a commonly owned, issued U.S. Patent, entitled EQUAL ORDER METHOD FOR FLUID FLOW SIMULATION, filed on even day herewith by the same inventive entity. The referenced application is incorporated herein by reference in its entirety.

[0002] This invention relates to a method for simulating the behavior of a fluid during a casting/molding process, in particular, to a numerical method for simulating incompressible, viscous melt flow in a die cavity.

[0003] In manufacturing casting/molding processes, for product quality of the molded article, it is of importance to optimize the mold design and the process conditions. Conventionally, such optimization is attained after many trial and error experiments which are repeated on the basis of the rule of thumb of an expert. With the advance in computer technology, it has become possible to analyze flow behavior of a material during a casting/molding process by computer simulation. When realistically modeled, the results from such simulation may be used to predict the casting behavior, thereby improving the efficiency of product and mold-design, and optimization of processing conditions.

[0004] Conventional methods for simulating fluid flow in a casting/molding process generally involve applying the equations of motion, continuity and energy to a discretized model of the molded part and performing numerical calculations thereof to obtain velocity, pressure, and temperature distributions along the processing line. In these simulations, the thermal physical properties of the material, the operational characteristic such as temperature of the material, temperature of the mold and injection speed, etc., are required input for the simulation. These inputs are generally assumptions made based on experimental know-how obtained from repeated comparisons between the analytical results and the actual moldings. The conventional models generally provide no means for determining whether or not input conditions are appropriate. It is therefore necessary to supply realistic boundary conditions before an accurate simulation can be achieved. Conventional fluid flow models lack many refinements which would enable better estimation of the boundary conditions, and hence often fail to achieve an accurate simulation.

[0005] One refinement which would be desirable, is incorporating the flow in the shot sleeve into the model. In a traditional die-casting process, molten material is delivered from the reactor to a die through a gated shot sleeve. Pressure is applied via a ram or plunger to eject the molten material from the shot sleeve. Due to discontinuity created by the moving ram and the resulting instability in numerical models, many prior art models start the simulation at the gate entrance and make the assumption that the pressure or velocity applied by the ram is uniform across the gate surface. In addition, prior art simulations assume the fluid contained within the shot sleeve has a uniform pressure and temperature. In practice, the molten metal in the shot sleeve always has some temperature variation, and the pressure applied to the billet across the ram face is non-uniform. More importantly, the thickness and temperature of biscuit remains in the shot sleeve after the cavity is filled has a major effect on the effectiveness of pressurization, which is required to suppress the formation of shrinkage porosity. This non-uniformity of temperature pressure and material flow in the short sleeve, if not corrected, ultimately compromises the accuracy of the mass flux computation of a flow model. Incorporating flow in the shot sleeve as part of the flow model would correct this deficiency.

[0006] Another refinement, which would be desirable, is to incorporate the heat exchange between the die and the heat-transfer fluid (HTF) into the model. A process parameter common to all metal die casting or plastic injection molding processes, including semi-solid forming, is the temperature control of the die. The temperature of the die is generally maintained by HTF contained in the heat transfer line which is embedded in the body of the die. Heated oil or water can be used to pre-heat the die to shorten the start-up time and the corresponding scrap. During continuous production, the temperature of the HTF can be adjusted to maintain the cyclic thermal balance in the die. Many attempts have been made to incorporate the heat transfer between the die and the HTF in a flow model. The most vigorous and resource-demanding approach is to use full 3D model to simulate the flow and convection heat transfer of the HTF inside the duct. In order to reduce the computational time, most of the existing models neglect the flow and temperature change of the HTF in the die. Instead, the user is required to select the element surfaces of the heat transfer lines and to assign, as thermal boundary conditions, a HTF temperature and a heat transfer coefficient. Since the HTF's viscosity changes with its temperature, which varies through the entire loop, it is very difficult for the user to determine the coefficient, as the Reynolds number and Prandtl number are not constant. A method which updates the heat transfer coefficient of the HTF as its temperature varies, such that the heat flux from the die to the HTF may be accurately determined, would be desirable.

[0007] A method for predicting the formation of shrinkage porosity would also be a desirable feature to be incorporated in a metal casting, e.g. semi-solid forming, simulation model. Shrinkage porosity formation is a common problem for metal casting process. When molten metal fills the cavity in a die which is at a lower temperature than the molten metal, the molten metal loses its thermal energy to the die and begins to solidify. As the liquid phase changes to solid, its density increases. If there is no additional material fed into the cavity, the net volume of the metal will become less than that of the cavity and shrinkage porosity will result. Depending on its size and distribution, shrinkage porosity could compromise the mechanical properties of a product, especially its elongation and fatigue strength, significantly. For pressure tight products, shrinkage porosity may randomly form multiple chained paths and cause leaking. Accurate predictions of shrinkage porosity are essential, when designing a structural component, e.g. control arm and knuckle, and high-pressure housing, e.g. air-conditioner housing and fuel-injection rail.

[0008] So far, the state-of-the-art casting simulation programs can only display potential locations of shrinkage porosity based on the hot spots (non-solidified metal surrounded by solidified metal) within the alloy's volume. In addition to temperature distribution, some of the programs also incorporate cooling rate and its empirical correlation to predict shrinkage porosity. However, the formation of shrinkage porosity is actually coupled not just with the thermal factors but also with the flow and pressurization. If the molten material in the cavity is in direct communication with a pressurized liquid metal source, e.g. the remaining liquid metal in the shot sleeve, additional material may be pushed into the cavity to compensate the volume shrinkage and to suppress the formation of shrinkage porosity. The required pressure depends on the viscosity of the solidifying metal, the alloy's temperature and the cavity's geometry. For accurate predictions and hence prevention of the formation of shrinkage porosity, a model that simulates both the thermal variation and the flow of the molten material would be required.

[0009] Another feature that is important to any casting or mold filling process is monitoring the cooling of the die after the part is ejected. Generally in a pressure casting process, after the metal solidifies in the cavity, the die opens to eject the part. A water-based lubricant is then sprayed on the mold walls for die release and to prevent drag mark or soldering. As the lubricant typically contains 98% water, it removes most of the thermal energy transferred to the die from the alloy. In order to predict the die temperature and its heat transfer with the alloy, the cooling effects of the lubricant spray in the model must be incorporated. However, due to the severe process conditions, there is no realistic method to determine the die temperature and hence the rate of heat removal by the lubricant. In addition, there are generally multiple nozzles spraying lubricants from different angles. Sometimes, these nozzles may move back and forth to spray lubricant more uniformly. Furthermore, depending on the type of nozzle, distance from the mold surface and spray pressure, the resulting distribution of cooling effects is not uniform. As a result, it is very difficult for user to define the corresponding cooling coefficient on the mold surface, especially one with complex shape. In some of the existing program, the user has to manually select and assign cooling coefficients on those elements surfaces that may be covered by the spray. In cases where that the nozzle moves, it would take too much time to set up the cooling coefficients. A model that computes the cooling efficiency from the geometry of the spray nozzle and its motion would be desirable.

[0010] A model to predict the location of a mend line and trapped air would also be a desirable feature for incorporation in a casting or molding simulation. When molten material fills a cavity, the melt front may split, e.g. to go around a core or a hole, and re-join. As the melt front is in contact with die steel and air, it may become too cold and/or contaminated by the residual lubricant or moisture. In addition, if the die is inappropriately designed, the cavity may not be filled progressively and air may be trapped inside the molten metal. As a result, wherever melt fronts meet, the bonding on the interface may be weakened due to solidification, trapped air or impurities.

[0011] Although the user can display the filling pattern from various angles to search for mend lines, it is very time consuming and not reliable, especially if the part is large and complicated. Furthermore, after the melt fronts meet, they could continue to move if the cavity was not completely filled. As a result, the final position of the weakened bonding interface may be different from the original mending location. In addition, the strength of bonding could be improved if there is a vigorous mixing of the material across the interface. It is very important in the design stage to be able to identify these potential casting defects quantitatively so that preventive design can be incorporated early in the design stage. Otherwise, the processing window would become very tight and it would be difficult to maintain consistent quality.

[0012] Therefore, there is a need for a new method which incorporates the aforementioned refinements such that incompressible, viscous fluid flow behavior in a casting/molding process can be accurately modeled by computer simulations.

[0013] Accordingly, an object of the present invention is to provide an improved method for simulating incompressible, viscous melt flow in a casting/molding process.

[0014] Related objects and advantages of the present invention will be apparent from the following description.

[0015] The present invention provides a method for simulating thermal flow during a molding/casting process, which satisfies the Navier-Stokes equations and the energy equation. This method incorporates five flow and thermal modules which simulate flow in shot-sleeve module, heat exchange in heat-transfer line, die cooling by lubricant, shrinkage porosity formation, and mend line formation in the casting/molding process. The incorporation of these modules provides better prediction of the boundary conditions and hence more accurate simulation. The modules may be incorporated in the overall casting model in any combination.

[0016] One aspect of the invention is a method for simulating thermal flows for casting/molding processes, comprising: (a) selecting simulation modules for incorporation in an overall simulation sequence, said modules including a Flow in Shot Sleeve Module, a Shrinkage Porosity Prediction Module, a Heat Transfer Fluid Line Module, and a Die Lubricant Cooling Module; discretizing flow domain into elements having nodes and fluxing surfaces, inputting initial parameter value and process conditions; (b) incrementing the time by an incrementing time step; (c) updating material balances; (d) forming and solving momentum equations, and updating velocity and pressure fields; (e) forming and solving energy equations, and updating temperature field; (f) checking if parameters have converged, returning to step (d) if they have not; (g) determining whether an end-simulation event has occurred, returning to step (c) if it has not; and (h) deciding whether to execute post-processor simulation modules, said post-processor simulation modules including a Mend Line Prediction Module.

[0027] For the purposes of promoting an understanding of the principles of the invention, reference will now be made to the embodiment illustrated in the drawings and specific language will be used to describe the same. It will nevertheless be understood that no limitation of the scope of the invention is thereby intended. Any alterations and further modifications in the illustrated device, and any further applications of the principles of the invention as illustrated therein being contemplated as would normally occur to one skilled in the art to which the invention relates are also included.

[0028] An embodiment of the present invention is a numerical method for computer simulation of the behavior of an incompressible, viscous fluid during a casting/molding process. The illustrated method is based on classical fluid dynamic models, and numerically (through the use of a computation device, generally a computer) obtains solutions to the governing equations at user defined time intervals. The illustrated method has general application to all thermal flows in casting/molding processes.

[0029] The method improves upon prior art methods by provides five novel modules: a Flow in Shot Sleeve Module, a Heat Transfer Line Module, a Die Lubricant Cooling Module, a Shrinkage Porosity Prediction Module, and a Mend Line Prediction Module, each of which simulates an aspect of the casting/molding process. Depending on the objective of simulation, engineer optionally selects to incorporate all or just some of the modules in a simulation. These modules provide realistic boundary conditions for solving the governing fluid dynamic equations. The sequential operations of each module is described in individual flow charts, separate from the overall simulation sequence. In these flow charts, the modules are depicted as complete sub-routines, mainly for illustrative purposes. However, it should be apparent to those skilled in the art that in an actual simulation where these modules are incorporated, the operations of the modules are integrated, in varying degrees, into the overall simulation sequence, and may become indistinguishable. The modules will be described in detail later.

[0030] The illustrated embodiment of the present invention employs a control volume-fixed mesh element model. The flow velocity and other physical parameters are calculated at the nodal points and may also be at interfacial surfaces. Although moving mesh methods are generally able to track the motions of the free surface more accurately than the fixed mesh methods, those approaches still have problems whenever interfaces come into each other or the die wall. A fixed mesh method is, therefore, preferred. However, a moving mesh model, which has overcome the moving-interfacial problems, may be used. In addition, any type of grid pattern in which the grid can be split into uniquely defined control volumes and that each of these control volumes can be associated with a unique node or grid point may be used. Preferably, the flow domain is discretized into a finite element mesh. The geometric flexibility of finite elements provides more accurate representation in fitting irregular domains and in providing local mesh refinement.

[0031] The illustrated embodiment of the method is based on a mesh composed of a three-node, two dimensional, triangular finite element similar to that shown in FIGS. 1 and 2. Referring now to FIG. 1, the triangle finite element 10 has three vertex nodes 12. Thermal physical parameters of the flow, e.g. temperature, velocity and pressures are evaluated at nodes 12 of element 10. Finite element 10 can be divided into three sub-control volumes (SCVs) 18 by joining the mid-side points 14 to the centroid 16. Each of sub-control volumes 18 is bounded by four fluxing surfaces 30 where material can flow from one sub-control volume to another. Typically, a number of sub-control volumes 18 surround a vertex node 12 and conveniently defines a control volume (CV) for that node. The sub-control volumes 18 that are attached to a particular node i, are collectively known as the control volume of node i. For example, as shown in FIG. 2, vertex node 22 is surrounded by five sub-control volumes 23–27 and the shaded area represents control volume 20 of vertex node 22. While it is illustrated that control volume 20 is defined by the sub-control volumes 23–27 attached to vertex node 22, control volumes can be composed of one or more sub-control volumes depending on the geometry of the flow domain. For a given mesh, a flexible combination of the SCVs is much better than solely relying on CVs in defining the irregular shape of the free surface of the flowing fluid.

[0032] For modeling a three-dimensional problem, where a three-dimensional mesh is required, a mesh composed of tetrahedral finite element is preferred. As shown in FIG. 3, tetrahedral element 40 has four vertex nodes 42 and four surfaces 43. Finite element 40 can be divided into four sub-control volumes (SCVs) 48 by joining the surface centroid 44 (centroid at each surface of the tetrahedral element 40) to the volume centroid 46 (the centroid of the tetrahedron). Each of the sub-control volumes 48 is a hexahedron having six fluxing surfaces.

[0033] While triangular and tetrahedral finite element are illustrated and preferred, it should be understood that other geometric shaped finite elements or a combination thereof may be used, and linear finite elements are preferred.

[0034] In this method the governing equations (1)–(3) can be expressed as:

• [0035] Continuity equation $\begin{array}{cc}\frac{\partial \rho }{\partial \tau }+\frac{\partial \left(\rho u\right)}{\partial x}+\frac{\partial \left(\rho v\right)}{\partial y}+\frac{\partial \left(\rho w\right)}{\partial z}=0& \left(1\right)\end{array}$
• [0036] Momentum Equation $\begin{array}{cc}\rho \frac{\partial V}{\partial \tau }+\rho \left(V·\nabla \right)V=-\nabla P+\rho g+\mu {\nabla }^{2}V& \left(2\right)\end{array}$
• [0037] Energy Equation $\begin{array}{cc}{\mathrm{pC}}_p\frac{\partial T}{\partial \tau }+\rho C_p\left(V·\nabla T\right)=\nabla \left(k\nabla T\right)+q^{″}+\Phi & \left(3\right)\end{array}$
  where
  • [0038] τ=time
  • [0039] ρ=density
  • [0040] P=pressure
  • [0041] V=the velocity vector
  • [0042] μ=viscosity
  • [0043] g=gravity components
  • [0044] T=temperature
  • [0045] C_p=heat capacity
  • [0046] Φ=viscous heating source term
  • [0047] q″=heat flux

[0048] Referring to FIG. 4 which shows the flowchart of the overall simulation sequence 100, according to an embodiment of the present invention.

[0049] The method starts at 1A where the mold is discretized. A finite element mesh of the flow domain may be created by any conventional procedures. For example, U.S. Pat. No. 5,010,501 issued on Apr. 24, 1991 to Arakawa and titled “Three-Dimensional Geometry Processing Method and Apparatus Therefor”, and U.S. Pat. No. 5,553,206 issued on Sep. 3, 1996 to Meshkat and titled “Method and System for Producing Mesh Representations of Objects”, both disclose suitable procedures.

[0050] At 1B, the matrix representing the mesh, the thermal physical properties of the fluid, initial values of the variables and other boundary conditions (e.g., temperature, pressure, density, velocity, the fill status of the SCVs, CVs, etc.) are input. The user also selects the specific modules to be incorporated in the simulation and the required boundary conditions necessary for the operation of these modules are also input.

[0051] At the beginning of the simulation, the master clock is initialized that t=0, 1C. At each subsequent time step, the time variable t, is incremented by a user defined incremental time interval Δt, 1D. The user may fix Δt or allows it to vary. It is preferred that Δt is variable. A variable Δt permits control of the incremental changes of the dependable variables in each time step, such that the number of iteration necessary to achieve convergence is even throughout the simulation.

[0052] At 1E, the status of each sub-control volume is determined: as filled, partially filled, or empty. The instantaneous volume of melt in a sub-control volume is the sum of its initial volume and the net-influx of melt from the beginning until the time step being analyzed. The mass flux is computed from the fluid density, velocity, the fluxing area and time. The amount of melt in a nodal control volume is the sum of the melt in all the attached sub-control volumes. Fill fraction volume is easily computed by using simple book keeping on each control volume. In the illustrated embodiment of the method, the mass flux contributions around a nodal control volume j is computed using a gauss point velocity field (u_g, v_g, w_g) that:
J_j=ρ(u_gω_xj+v_gω_yj+w_gω_zj)  (4)
where J_j is the mass flux around a nodal control volume j, ρ is density, and ω is the fluxing surface area.

[0053] The development of the gauss point velocity field is disclosed in an application entitled “Equal Order Method for Fluid Flow Simulation”, filed on even day with the subject application by the same inventors, and hereby incorporated by reference in its entirety. Computing mass flux using the gauss point velocity field conserve mass to the truncation limits of the computation device (machine round-off levels). In addition, the method also uses a weighting scheme to compute mass flux to ensure global mass balance is conserved. For mass fluxes that are filling a nodal control volume, the upstream fill fraction is used to weight the filling flux. If a mass flux is draining a nodal control volume, then the downstream fill fraction is used to weight the flux. This weighting scheme is much like the donor-acceptor algorithm used by SOLAVOF.

[0054] While using the gauss point velocity field to compute mass flux is preferred, other scheme of computing mass flux, e.g. the SIMPLE and SIMPLER methods, may be used with the present method.

[0055] At this stage, the position of the free surface is also located. The volume fraction of the melt, f_m, is used to infer the position of the free surface. f_m, is defined as:
f_m=V_m/V  (5)
where

• • [0056] V and V_m denote the volume of the control volume and the volume occupied by the melt, respectively.

[0057] In an element where the melt front exists, the volume fraction of the melt, f_m, on some of the nodes would be 1 and others less than 1. The free surface is defined as the surface where the volume fraction of melt, as interpolated from the nodes, is 0.5. However, since the exact position of the free surface is not known and some mass errors can accumulate. To avoid this error, the mass in the computational domain is re-normalized after each time step to ensure exact mass conservation.

[0058] After the filling status of the sub-control volumes are determined, the trapped air pockets are identified as a single or a group of empty sub-control volumes surrounded by filled or half filled sub-control volumes. These trapped air pockets are named and are assigned atmospheric pressure as their initial pressure. The position and melt fraction of these trapped air pockets are tracked and updated with each time step. The pressure within the trapped air pocket is calculated from the ideal gas law using the total mass flux into the pocket during a specific time interval. The updated pressure becomes the boundary condition for any in-filling event at the next iteration.

[0059] A computational domain is defined after updating the filling status of the nodal control volumes and the melt front. The purpose of defining a computational domain is to conserve computation time, only control volumes within the computational domain will be considered in the calculation. In the illustrated embodiment, the border of the computation domain is defined by joining nodes which are adjoining free-surface nodes, which have at least one filled sub-control attached. However, other criteria of defining a computational domain may be used.

[0060] At IF and 1G, the method determines whether the user optionally selects to incorporate the Flow in Shot Sleeve Module 200 and the Shrinkage Porosity Prediction Module 300 in the overall simulation. If the user has chosen to incorporate these modules, the step sequence of these modules, 200 and 300 will be integrated in the overall solution sequence, such that the momentum equations will be formed and solved with the appropriate terms and boundary conditions. The step sequences of modules 200 and 300 will be discussed in detail later.

[0061] At 1H, the momentum equations and the pressure equations are formed and solved. These equations govern the conservation of mass and momentum. Each momentum equation is made up of terms defined by the user. The terms commonly included for viscous, incompressible flow are inertia, advection, diffusion, body force and a pressure gradient. Using CVFEM techniques, an elemental matrix can be computed to approximate each term. These matrices are then assembled into a global momentum matrix and solved iteratively. In the illustrated embodiment, the global momentum equations are formed and solved by the Equal Order method.

[0062] The Equal Order method is a numerical method for solving the coupled momentum and mass conservation equation. The method defines a gauss point velocity vector as the average of the nodal velocity vectors. In essence, the expression of gauss point velocity is in the same order as the pressure gradient over each of the elements. Because the gauss point fluxing vector field is “centered” inside the pressure field defined at the nodes, no pressure checker boarding will occur. Additionally, because the gauss point velocities and pressure fields are evaluated from the same number of nodal positions, the accuracy of the solution is preserved. For linear triangular elements and linear tetrahedral elements, mass can be conserved to machine round-off levels. Using a two dimensional flow as illustration, the assembled momentum equations are written as: $\begin{array}{cc}\sum _{j=1}^N{a}_{\mathrm{ij}}{u}_j=-{\int }_A\frac{\partial P}{\partial x}dA\mathrm{and}& \left(6\right)\\ \sum _{j=1}^N{a}_{\mathrm{ij}}{v}_j=-{\int }_A\frac{\partial P}{\partial y}dA& \left(7\right)\end{array}$
where the a_ij terms make up the influence coefficient matrix. Pulling the pressure gradients outside the area integrals, these equations are recast as: $\begin{array}{cc}{u}_i={\hat{u}}_i-K_i\frac{\partial P}{\partial x}& \left(8\right)\end{array}$ $\begin{array}{cc}{v}_i={\hat{v}}_i-K_i\frac{\partial P}{\partial y}& \left(9\right)\end{array}$
where, $\begin{array}{cc}{K}_i=\frac{1}{a_{\mathrm{ii}}}\int dA,& \left(10\right)\\ \hat{u}=-\frac{1}{a_{\mathrm{ii}}}\sum _{i\ne j}{a}_{\mathrm{ij}}{u}_j& \left(11\right)\\ \hat{v}=-\frac{1}{a_{\mathrm{ii}}}\sum _{i\ne j}{a}_{\mathrm{ij}}{v}_j& \left(12\right)\end{array}$
and the nodal hat velocities, û_i and {circumflex over (v)}_i, are composed of nodal velocities u_j and v_j of neighboring nodes and contain no pressure term. The gauss point hat velocity û_g, {circumflex over (v)}_g and the gauss point coefficient term K_g are defined as: $\begin{array}{cc}{\hat{u}}_g=\frac{1}{N_e}\sum _{j=1}^{N_e}{\hat{u}}_j,& \left(13\right)\\ {\hat{v}}_g=\frac{1}{N_e}\sum _{j=1}^{N_e}{\hat{v}}_j,& \left(14\right)\\ K_g=\frac{1}{N_e}\sum _{j=1}^{N_e}{K}_j& \left(15\right)\end{array}$
and the gauss point velocities are expressed as: $\begin{array}{cc}{u}_g={\hat{u}}_g-K_g\frac{\partial P}{\partial x}& \left(16\right)\\ v_g={\hat{v}}_g-K_g\frac{\partial P}{\partial y}& \left(17\right)\end{array}$
In this expression for the gauss point fluxing velocities, u_g and v_g, the pressure gradients, ∂P/∂X, ∂P/∂Y, are constant over the fluxing surface. For consistency, the terms u_g, v_g and K_g must all be constant over the fluxing surface as well.

[0063] By inserting the expression for the gauss point fluxing velocities, u_g and v_g, into the continuity equation (1), a Poisson equation for pressure is obtained. The pressure equation is expressed as: $\begin{array}{cc}\nabla k\nabla p=\frac{\partial \hat{u}}{\partial x}+\frac{\partial \hat{v}}{\partial y}& \left(18\right)\end{array}$
Integrating this equation over the area of the fluxing surface of the eelement, the pressure field can be expressed as: $\begin{array}{cc}{\int }_A\nabla k\nabla \mathrm{PdA}={\int }_A\left[\frac{\partial \hat{u}}{\partial x}+\frac{\partial \hat{v}}{\partial y}\right]dA& \left(19\right)\end{array}$
Using Green's theorem, these area integrals are converted to line integrals, $\begin{array}{cc}{\int }_{\sigma }\left[k_g\frac{\partial P}{\partial x}\right]dy+\left[k_g\frac{\partial P}{\partial x}\right]dx={\int }_{\sigma }{\hat{u}}_{g}dy-{\hat{v}}_{g}dx& \left(20\right)\end{array}$

[0064] Using the gauss point hat velocities û_g, {circumflex over (v)}_g and K_g previously defined, the Poisson pressure equation can be assembled and solved to determine the pressure field.

[0065] The nodal velocities, u_i and v_i and the gauss point velocities, u_g and v_g, are computed by inserting the newly derived pressure field and û_i, {circumflex over (v)}_i, K_i, û_g, {circumflex over (v)}_g and K_g in the momentum equations.

[0066] While the Equal Order method is preferred, it should be understood that other methods which derive and solve the discretized form of the momentum equation may be used without deviating from the spirit and scope of this invention. Examples of the other known methodology are the SIMPLE and SIMPLER methods of Baliga and Patankar, (in “A Control Volume-Finite Element Method for Two Dimensional Fluid Flows and Heat Transfer”, in Numerical Heat Transfer, 6, 245–261 (1980)), the Galerkin's method, and the Newton Raphson method.

[0067] At 1I, the method determines whether the user has elected to incorporate the HTF Line Module 400 in the simulation. If the user made the election, the step sequence of module 400 will be inserted before 1J of the overall simulation sequence 100. The sequential steps of the HTF Line Module 400 will be discussed later. The heat flux from the die to the heat transfer fluid, a result of module 400, will be used to predict the die temperature.

[0068] After solving the pressure and velocity field, the method proceeds to solve the temperature fields for both the die and the melt by solving the energy equation (3). The energy equations are made up of terms that are defined by the user. Preferably, the energy equation for the melt includes terms for inertia, advection, diffusion, volumetric heat source, viscous heating, and coincident interface heating. Additional, Dirichlet and Neumann boundary conditions should be included as well. The energy equation for the die preferably includes terms for conduction and radiation. Using CVFEM techniques, an elemental matrix is computed to approximate each term. These matrixes are then assembled into a global energy matrix. After the assembly process, the left hand side can be represented by a coefficient matrix operating on a temperature solution vector. The right hand side is made up of a source term vector. This set of discretized form of the energy equations is solved to yield the temperature at each node. Fraction of solid information is computed using a temperature versus fraction of solid curve. The curve can be empiracally determined.

[0069] At 1K, the method checks for convergence of the dependent variables. An iterative process is said to have converged when further iterations will not improve the accuracy of the dependent variables. In practice, the iterative process is terminated when the desired accuracy is obtained. Although the computed velocity field at the gauss points conserves mass to machine round-off at every iteration, the mometum equations are non-linear and therefore iterations are required to refine the momentum balance. In this technique, the flow equations and energy equations are iterative solved until all velocity and temperature values converge to a user supplied tolerance. If either of the nodal or gauss point velocities or temperatures, Ti has not converged, the method returns to 1H where the momentum and energy equations are solved and iterated again using the updated nodal and gauss point velocities and Ti. If the velocity and temperature fields have converged, the method then proceeds to the next stage.

[0070] At 1L, the decision is then made whether or not to store the calculated parameters and which parameters to store. These decisions are based on the requirements of the user of the method—results can be stored every time step, every other time step, or some other increment desired by the user. Preferably, the parameters being stored include: time, velocities, pressure, die and melt temperatures, solid fraction, shrinkage porosity, heat transfer fluid temperature and the free surface position. These values may also be output.

[0071] At 1M, the method checks for control trips defined by user. If a trip condition has been reached, the method will follow a user-defined sequence of events. In most cases, trip conditions are used to signal the computer program to make boundary condition changes or to begin or end the use of a particular module.

[0072] At 1N, the method performs an energy audit, such that the heat flux of all the components may be integrated over the time step. The energy audit provides a quantitative understanding of the net energy exchange between components such that the user may assess whether the appropriate process settings were provided for the simulation.

[0073] At 1O, the method checks for end of cycle conditions. If the simulated time has reached the end of cycle time, then the initial conditions are re-applied to the model and a new cycle is simulated. This allows for die heat-up transients to be studied.

[0074] At 1P, the method allows the user to define the next incremental time Δt. If user does not define a incremental time Δt, the method will calculate the next time interval based on the number of iterations performed to obtain convergence for the present time step.

[0075] At 1Q, the method then checks to see if the stop flow simulation time has been reached. The stop flow simulation time is a user-defined value. If the stop time has not been reached, the method returns to incrementing the time at 1D, and proceeds with another cycle. If the stop flow simulation time has been reached, then the mold filling portion of the simulation is completed. However, the thermal portion of the simulation continues if the user elected to incorporate the Die Lubricant Cooling Module 500 in the simulation at 1R, the operation sequence of module 500 will be incorporated into the overall simulation sequence 100. The sequential operations of the Die Lubricant Cooling Module 500 will be discussed later.

[0076] At 1S, the user has the option to simulate for mend line formation. Mend line prediction is a post-processing activity and does not interact with the computational cycle of the overall simulation sequence 100. The sequential operations of the Mend Line Prediction Module 600 will be discussed later.

[0077] Simulating the flow in shot-sleeve model established the initial pressure, temperature and velocity vector field of the flow before the mold-filling event occurs. The pressure field can be used to establish the dwell time, dwell pressure and the thickness of the biscuit (material left in shot sleeve after the mold is filled) required for maintaining the continuity of mass flux. Dwell time and dwell pressure are important process parameters for avoiding the formation of shrinkage porosity. Furthermore, the module allows the user to set initial temperature field of the billet and monitor the movement of hot and cold spots as the mold cavity is filled.

[0078] When the program runs, the position of the ram in each time interval is monitored. The ram displacement is used to determine the amount of volume of molten metal displaced from each control volume adjacent to the ram. This volume of metal is then incorporated in the continuity and momentum equations to calculate the pressure and velocity fields. The temperature field is determined from the energy equation.

[0079] Referring now to FIG. 5 which shows the flowchart of an embodiment of the Flow in Shot Sleeve Module 200 of the present invention. The sequential operation of Flow in Shot Sleeve Module 200 is highly integrated into the overall simulation sequence 100. While the module 200 is depicted in FIG. 5 as a complete, independent routine, it is so depicted for conceptual understanding and demonstrative clarity. It is contemplated and understood that it will be readily apparent to those skilled in the art that the sequential operation of Flow in Shot Sleeve module 200 is interspersed between, and relies upon the operation of the overall simulation sequence to function.

[0080] The module starts at 2A with the inputting of the geometric model of the shot sleeve and the ram. The shot sleeve and the ram are then meshed into a plurality of finite elements and the matrix representing the mesh is entered. At 2B, the user may define a set of process parameters and boundary conditions comprising: the initial position of the ram, the velocity profile of the ram, the volume fraction of the melt in the mold reaching when the speed control switches from velocity to dwell pressure, the dwell pressure profile of ram, the length of time that the ram is to be held at dwell pressure, initial volume and temperature distribution of billet. The thermal physical properties of the material are also input.

[0081] At 2C, the master clock advances by an incremental time Δt. The time clock is the same as 1D of overall solution sequence 100, and the incremental time Δt is also the same as the overall solution sequence 100.

[0082] At each time step, the position of the ram is updated, 2D. As the ram always moves in its axial direction, the ram position is determined from the user define ram velocity profile and the time increment Δt.

[0083] From knowing the ram position and the stored matrix of the shot sleeve, the method identifies the control volumes which have been crossed by the incremental advancement of the ram using simple book keeping techniques (2E). The control volumes of the shot sleeve, which are occupied by the ram, are deleted from the computation domain.

[0084] Updating pressure and velocity fields at 2F is the same operation as 1H of the overall simulation sequence 100. The incorporation of Flow in Shot Sleeve Module 200 modifies the boundary conditions for the computation of the velocity and pressure fields. The volume of molten metal displaced by the advancing ram is computed from the geometry of the shot sleeve and the incremental ram positions. This volume of molten metal is equated to the mass flux across the fluxing faces of the control volumes. The mass flux is inserted into the continuity equation to form a Poisson equation for pressure. The pressure equation is solved to obtain the pressure field.

[0085] Using the updated pressure field, the global momentum equations are solved for the velocity field. The momentum equations are formed and solved by the Equal Order method as previous described in 1H of overall solution sequence 100. It should be understood that while the Equal Order method is preferred, other method of forming and solving the momentum equation might be used without deviating from the spirit and scope of the invention.

[0086] Updating temperature field at 2G is the same operation as 1J of the overall simulation sequence 100. Again, the incorporation of Module 200 modifies the boundary conditions for the computation of the temperature field. In regions where coincident elements remain at their initial positions, aligning-coincident heat transfer is assumed, which is similar to that described in 1J of overall solution sequence 100. Non-coincident (or non-aligning) heat transfer is also present and results from the motion of the ram such that the elements of the shot sleeve or the melt in contact with the ram are constantly changing. The method accounts for this unique heat transfer by vigorous book keeping of the interfacial boundary.

[0087] The method checks for convergence of the velocity and temperature fields at 2H (same as 1K of overall simulation sequence 100). If the velocity field or temperature field has not converged, the method returns to 2F for another iteration using updated pressure, velocity and temperature values. If the velocity and temperature has converged, the method checks whether the criteria for the ram control to switch from speed to dwell have been reached, 2I. Preferably, the criterion for switching is the fill fraction of the mold. The user may select the value of the criterion; preferable, it is set at 99% fill. Alternatively, the user can also choose a time or ram position to switch ram speed control into pressure control.

[0088] When in dwell, the pressure of the ram forms the boundary condition of the elements immediately adjacent to the ram, 2J. The pressure field is updated. The motion of the ram is now updated based on the amount of material pushed out from those adjacent control volumes. The method then checks whether end time for dwell has been reached, 2K. The end time for dwell is selected by the user. If end time for dwell has not been reached, the method returns to 2C and advances a time increment Δt and the pressure, velocity field are updated. At dwell, the driving force for melt flow is pressure rather mass flux. If end time for dwell has been reached, the method continues to 1R of overall solution sequence 100 (2L).

[0089] The porosity prediction module 300 uses a coupled thermal-flow model to predict the formation of shrinkage porosity. First, the thermal-flow model provides a prediction of temperature, pressure and velocities within the region occupied by the metal and the die. Under the constraint of the conservation of mass, a region where alloy is becoming denser over time will cause a very low pressures in the momentum equation if the flow resistance in the surrounding matrix is large enough to retard material from filling in. As the local pressure approaches vacuum, the region would not be able to suck any additional material from the surrounding matrix to compensate for its shrinkage. Any further shrinkage will then produce porosity in that region. Therefore, shrinkage porosity can be calculated directly by incorporating porosity and density variation in the mass conservation equation. This process can be undertaken at each time level to produce a time history of porosity at each position in the model.

[0090] Referring now to FIG. 6 which shows the flowchart of an embodiment of the shrinkage porosity prediction module 300 of the present invention.

[0091] At 3A, the momentum equations are discretized and expressions for the velocity components are formed as described in 1H of overall simulation sequence 100. The continuity constraint is modified to include the effects of porosity. The modified continuity constraint is expressed as: $\begin{array}{cc}\frac{\partial \beta u}{\partial x}+\frac{\partial \beta v}{\partial y}+\frac{\partial \beta w}{\partial z}+\frac{\partial \beta }{\partial t}=0& \left(21\right)\end{array}$
where

• • [0092] β=(1−φ)ρ apparent density of the alloy
  • [0093] x, y, z=three components of the coordinate system
  • [0094] t=time
  • [0095] u, v, w=three velocity components
  • [0096] φ=volumetric fraction of porosity (0— no porosity and 1–100% porosity)
  • [0097] ρ=theoretical density of the alloy as function of temperature

[0098] The velocity components are inserted into the modified continuity constraint to form the Poisson equation for pressure. This equation is solved to determine the pressure at each node within the computational domain.

[0099] At 3B, the pressure field is analyzed. Shrinkage porosity forms or increases when the control volume's pressure decreases to vacuum. The control volumes having negative pressures are identified and the pressure of those control volumes are re-set to zero, 3C. Using the reduced pressure gradient, the momentum equation is solved and the velocity field updated, 3D. At 3E, the volumetric net inflow rate, {dot over (v)}, in each control volume with negative pressure is evaluated. The reduced pressure gradient yields a smaller net flow rate than when the pressure in the specific control volume was negative. The smaller net flow rate may be rationalized as accounting for the space taken up by the voids.

[0100] At 3F, the increment in porosity within the control volume is computed by taking the difference between the densification rate and the smaller net flow. The following equation is used to estimate the increase of porosity within the control volume: $\begin{array}{cc}d\varphi =\frac{1}{{\rho }_s}\left(\frac{\partial \rho }{\partial \tau }-\frac{\stackrel{.}{v}}{V}\right)d\tau & \left(22\right)\end{array}$
where

• • [0101] ρ_s is the theoretical density at solidus temperature
  • [0102] ρ is the theoretical density
  • [0103] V is the volume of the control volume
  • [0104] {dot over (v)} is the net flow rate of material coming into the control volume
  • [0105] dτ is the incremental time step.

[0106] At 3G, the apparent density β, where β=(1−φ)ρ, at each node is calculated using the updated porosity value φ. The updated apparent density β is used in computing the velocity and pressure fields in the next time step. A time history of shrinkage porosity at each control may therefore be established. At 3H, module 300 continues to 1I of overall simulation sequence 100.

[0107] The HTF line module 400 simulates the heat flux between the die surface and the heat transfer line, such that the temperature field may be determined accurately and processing conditions may be set properly. FIG. 7 shows a schematic diagram of the heat transfer line 60 in relationship to the die elements 61. The entire heat transfer line 60 is divided into a number of linear segments k_j, where j=1,2, . . . N. Each segment k_j is associated with a known set of die element surfaces 62 at where the heat transfer from the die to the cooling channel is accounted. Each segment k_j is bounded by an upstream node l_j-1 and a downstream node l_j, named in relation to the flow direction m. The average temperature distribution within each segment k_j is linearly interpolated from the nodal temperature at its two ends. As the HTF lines are typically drilled with circular cross section, the axisymmetric velocity and temperature distribution of the HTF can be calculated based on the well-established analytic solutions. This model, thus, simplified a three-dimensional problem to a one-dimensional problem in which the solution is manageable.

[0108] The governing equation for the this heat transfer is an implicit form of a one-dimensional energy equation: $\begin{array}{cc}\frac{A_{\mathrm{flow}}L\rho C_p}{\mathrm{\Delta \tau }}\left[\frac{T_j+T_{j-1}}{2}-\frac{T_j^0+T_{j-1}^0}{2}\right]=-{\mathrm{mC}}_p\left(T_j-T_{j-1}\right)+\sum _{k=1}^{N_f}{h}_k{A}_k\left(T_{\mathrm{mk}}-\frac{T_j+T_{j-1}}{2}\right)& \left(23\right)\end{array}$
where

• • [0109] A_flow=the cross sectional area of the HTF line
  • [0110] Δτ=time step
  • [0111] h_k=convection heat transfer coefficient between the HTF and the die in segment k
  • [0112] A_k=circular area of the HTF line in segment k
  • [0113] T_mk=mean die temperature in segment k
  • [0114] T_j=temperature of the HTF at node j at the end of a time step
  • [0115] T^o_j=temperature of the HTF at node j at the beginning of a time step
    By defining the following terms, $\begin{array}{cc}I=\frac{\Delta \tau }{A_{\mathrm{flow}}L\rho C_p}& \left(24\right)\\ E=\sum _{k=1}^{N_j}{h}_k{A}_k,& \left(25\right)\\ E_T=\sum _{k=1}^{N_f}{h}_k{A}_k{T}_{\mathrm{mk}}& \left(26\right)\end{array}$
    an update of the cooling channel temperature, T_j, can be expressed as,
    T_j=(cbT_j-1)/α  (27)
    where,
    α=1+I(E+2{dot over (m)}C_p)  (28)
    b=1+I(E2{dot over (m)}C_p)  (29)
    c=2IE_T+T⁰_j+T⁰_j-1  (20)

[0116] The method of the present invention solve the above energy equation numerically using a tri-diagonal solver. At each time step, the Reynolds number (Re), Prandtl number (Pr), Nusselt number (Nu), and the corresponding heat transfer coefficient h_k on each element surfaces of the heat transfer line is calculated based on the user defined volume flow rate, thermal physical properties of the HTF, and the updated HTF temperature. The heat flux exchanged between the HTF and the die is calculated on each die element surface 62 from the heat transfer coefficient h_k, the local temperature difference between the HTF and the die, the element surface area and the time step. With the heat flux of the die element surfaces 62 within each line segment k_j, the HTF's temperature at the upstream node l_j-1, the volume flow rate and thermal physical properties of the HTF, the HTF's thermal energy, and setting the boundary condition that the temperature of the first node T₀=T_inlet, temperature at the downstream node T₁ can be calculated. By this approach, the analytical/empirical correlation between Re, Pr and Nu and the variation of the HTF's flow rate, temperature and viscosity can be incorporated to predict the heat transfer accurately without heavy increase in the computation steps.

[0117] Referring now to FIG. 8 which shows the flowchart of an embodiment of the Heat Transfer Line Module 400 of the present invention.

[0118] At 4A and 4B, the method identifies all of the die element surfaces 62 associated with the HTF lines from the imported finite-element mesh of the entire die and user specified element edges of the inlet and outlet. Then, these HTF element surfaces are grouped sequentially, based on their flow distance from the inlet, into a number of segment k_j where j=1, N. The number of segments, N, is determined by the user based on the total length of the HTF line.

[0119] At 4C, based on the initial temperature T_k of the heat transfer fluid at the beginning of a time step, the method updates the viscosity, μ, of the HTF at segment k_j, from a look-up table which contains empirically determined viscosity and temperature relationship for the fluid. From viscosity, Reynolds number, Re=ρvd/μ, and Prandtl number, Pr=C_pμ/k, are calculated. The Nusselt number, Nu, being a function of the Reynolds and Prandtl numbers can also be updated. The relationship between the between Re, Pr and Nu numbers can be analytically/empirically determined for a specific fluid.

[0120] At 4D, the heat transfer coefficient h_k in segment k_j is calculated from the Nusselt number using the following relationship:
h_k=Nuk/d
where h_k=heat transfer coefficient of the HTF in segment k

• • [0121] k=thermal conductivity
  • [0122] d=diameter of the line

[0123] At 4E, the value of T_mk is input into the calculation. T_mk is the average nodal temperature of the set of die elements 61 associated with segment k, of cooling line. The nodal temperatures of the die elements 61 are read from the temperature field of a previous iteration.

[0124] At 4F, the transport equation is built, a tri-diagonal matrix is formed for the HTF temperature. Given the boundary condition that T_j=0=T_inlet, the average temperatures at every node l_j is computed.

[0125] At 4G, the updated temperatures are then used as boundary conditions to calculate the heat flux between the die and the HTF to solve the temperature field in the cavity and the die at 1J of overall simulation sequence 100 at 4H.

[0126] The Die Lubricant Cooling Module 500 calculates the amount of heat removed from the die by spraying of lubricant after the part was ejected. In this module, the element surfaces covered by the lubricant spray are identified mathematically by emanating a spray cone from the nozzle's position toward the mold surface. By comparing the relative position of the element surfaces on the mold cavity with the spray cones, those element surfaces covered by the spray can be identified. If the spray nozzle moves, the position of the spray cone and the corresponding sprayed element surface will be updated for every time step. The cooling coefficient on an element surface is determined empirically. In actual set up, the user only needs to define the location, angle and motion path of each spray nozzle, nozzle type, spray pressure and flow rate. The program will automatically identify and assign a time-dependent convection-cooling coefficient at every time step on those element surfaces intercepted by the spray cone. With this function, the user can quickly change the spray nozzle, spray pattern and spray parameters to evaluate their impact on die temperature and product's quality.

[0127] Referring now to FIG. 9 which shows the flowchart of an embodiment of the Die Lubricant Cooling Module 500 of the present invention. After stop flow time for the flow portion of the simulation has been reached, the die opens and the part is ejected. The matrix of the open die is input for the Die Lubricant Cooling Module 500. At 5A, a local coordinate system is established using the nozzle location as the origin, the major and minor axis of the nozzle as the x and y coordinate and the axis extend from the nozzle to the die face as the z-axis. A transformation matrix is used to place each element of the three-dimensional mesh that makes up the mold wall surface into the local spray coordinate system. Alternatively, the cone can also be represented in the global coordinate of the imported finite-element mesh.

[0128] At 5B, the user is required to specify the physical properties and position of each spray nozzle, e.g. nozzle shape, the major and minor axis of the nozzle and its distance from the die surface, and the path of travel. In addition, the user has to specify the flow rate and the spray pressure. The user can vary the parameters to achieve different spraying conditions. For example, to change nozzle shape to accommodate a different spray pattern, or the spraying time. Preferably the user would set the time of the spray and travel path of the nozzle so that the lubricant would spray the entire mold surface.

[0129] The Die Lubricant Cooling Module 500 is integrated into the global energy balance of the simulation. The master clock at 1C of overall simulation sequence 100 controls the operation of this module, 5C.

[0130] At 5D, the nozzle's location in the local coordinate system is input. At 5E, from the shape of the nozzle, the major and minor axis of the nozzle, and the distance from the surface, the method simulates the shape of the coverage of the spray or the spray cone. It should be understood that the shape of the spray cone depends on the shape of the nozzle, which can be any shape. Preferably, the shape of the nozzle is circular or elliptical. At 5F, the elements facing the nozzle are identified. The identification is accomplished by a visual test, which the nozzle is placed at the eye position and viewed along the Z-axis of the local coordinate system towards the die. The die element faces which are in front of the spray and having a normal vector such that they face the spray are selected. Upon the element surfaces that are selected, a second test is performed to determine which element faces are actually exposed to the spray (i.e., within the cone of influence), instead of being hidden behind other die features, e.g., a moving core. The test is: $\begin{array}{cc}\mathrm{if}\frac{x_l^2}{a^2}+\frac{y_l^2}{b^2}<z_l^2,& \left(31\right)\end{array}$
then, the element face is determined to be inside cone of influence. The subscript, l, is used to denote the local coordinate system as defined by the spray.

[0131] At 5H, if a face is in the cone of influence, the face is projected through a hemispherical surface to compute the fraction of the total view space that this face occupies. The hemispherical projection technique is well-documented in radiation heat transfer literature (e.g., Heat Transfer—A Basic Approach, M. N. Ozisik McGraw Hill, 1985.) Once the projected area of an element face is known, the heat transfer coefficient at 5I may be look-up from a table and the heat flux from the element face may then be computed at 5J. The heat flux is inserted into the energy equation to update die temperature for the subsequent process.

[0132] The look-up table is constructed empirically. The following is an example of the construction of such a table. First the volume of spray reaching the projected area of the element face is computed by the expression:
Q_f=Q_sη_sf_hp  (32)

• • [0133] Q_f=volume of spray reaching mold face,
  • [0134] Q_s=Volumetric spray flow
    Where,
  • [0135] η_s=Spray carry out efficiency
  • [0136] f_hp=fraction of the view space

[0137] Then spray experiments can be conducted where a known volume of lubricant is sprayed on a heated flat plate. All thermal physical parameters of the flat plates are known. Preferably, the plate is of similar material as the die. The temperature drop resulted from the spray is measured, and the heat transfer coefficient computed. The volume of spray used versus heat transfer coefficients are stored. The table can then be used to look up heat transfer coefficient when the volume of spray reaching the element faces of the die is known. While a suggestion is made as how to establish a look-up table for heat transfer coefficient, other methods, empirical or numerical, of determining a heat transfer coefficient from the projected area of a die element face may be used.

[0138] After the heat flux and die temperature are calculated, the method checks whether a user defined stop spray time has been reached at 5K. If stop spray time has not been reached, the sequence returns to 5C, and time is incremented and the spraying sequence repeated. If stop spray time has been reached, the Die Lubricant Cooling Module 500 is completed. At 5L the simulation continues 1S of overall solution sequence 100.

[0139] The Mend Line Prediction Module 600 simulates the location of the mend line in the flow. Module 600 predicts the position of the mend from the data collected during a mold filling simulation. The module is executed after the completion of the simulation. It is a post-processing activity and does not interact with the computational cycle of the overall simulation sequence 100.

[0140] Along the free surface of the melt front, mathematical points are created to represent material particles. The position of each particle is tracked. At each time step, particles are advected to their new position based on the calculated melt velocity on the previous position. After the simulation is completed, the user can display the distribution of these particles which were assigned different colors or symbols denoting the time of birth (when they were added) or local temperature. In this way, mend lines can be identified.

[0141] Referring now to FIG. 10 which shows the flowchart of an embodiment of the Mend Line Prediction Module 600 of the present invention. At start, a clock is initialized, user sets start time and end time for the simulation. User may choose to predict the formation of the mend line only at a specific time period of the overall simulation at 6A.

[0142] The user may set the incremental time Δt between time steps. User may follow the same incremental time, Δt, used in the simulation where the data is collected. Alternatively, the user may define the time increment according to the resolution of the motion of the melt desired, 6B.

[0143] At 6C, the position and shape of the melt front in each partially-filled control volume is mathematically identified based on the volume fraction of the melt, f_m. The volume fraction of the melt is defined as: $\begin{array}{cc}{f}_m=\frac{V_m}{V}& \left(35\right)\end{array}$
where V_m is the volume of the control volume occupied by the melt

• • [0144] V is the total volume of the control volume
    In one embodiment of the present invention, the free surface of the melt front is defined at f_m=0.5. The elements where the melt front exists are first identified from the volume fraction of their nodes. In an element where the melt front exists, the volume fraction of the melt on some of the nodes would be larger than 0.5 and others less than 0.5. The volume fraction of melt in each control volume is updated based on the incremental net inflow of melt during the iterations in every time step.

[0145] At 6D, the module then determines whether the number of particles at the melt front is sufficient. The module determines the distance between the attached particles and the melt front, if the distance between a particle from the free surface of the melt front is larger than a user defined value, the particle is considered detached. If the number of particles on a melt front surface in an element is lower than the user-specified number density, the module will add new particles to the free surface. In addition, the method assigns different colors or symbols to track these newly added particles. If the number of particle at the free surface is judged sufficient, the module continues to 6G.

[0146] While a distance criteria is used for adding particles, other criteria, such as a range of fraction of liquid, particle density, etc., may be used to determine whether a particle is still attached to the melt front without deviating from spirit and scope of the invention.

[0147] At 6E, new particles are created/added as mathematical points along the free surface of the melt front. Each particle is assigned a position, x_p, y_p, z_p a temperature t_p^o and a velocity V_o^o which are equal to the velocity and temperature of the associated melt front, as interpolated inside an element from the nodes. Once created, the position, temperature and velocity of each particle are stored and will be tracked continuously, 6F.

[0148] At 6G, the position of all particles are updated. At each time step, the particle is advected to a new position. The new position is updated using the local velocity vectors.
x′_p=x_p+v_xΔτ
y′_p=y_p+v_yΔτ
z′_p=z_p+v_zΔτ
where (x′_p, y′_p, z′_p) is the latest position of the particle, and v_x, v_y, v_z are the velocity component of the particle. During this advection process, the particles are constrained such that they can not pierce the boundaries of the model. The new position is stored.

[0149] At each time step, the average position of the free surface within each element is updated. At each element that contains a free surface, the intersection of the f_m=0.5 surface with the edges of the element is computed. The resulting collection of points is used to define a polygon whose centroid is denoted by x_fs, y_fs, z_fs, 6H.

[0150] At 6I, the method determines the attach/detract status of each particle. The method compares the latest position of each particle (x′_p, y′_p, z′_p) to the average position of the free surface (x_fs, y_fs, z_fs). If the particle is within a user-defined tolerance of the free surface points, then the particle is defined as an attached particle. If the particle is outside the defined tolerance of the free surface points, then it becomes a detached particle. The detached particles are labeled with different color/symbol for easy visual identification.

[0151] At 6J, the particle locations are rendered to the computer's screen for visual analysis. The attached and detached particles are distinguished by color/symbols. Based on the display, the user can easily identify the mend line, where particles of different colors meet, its final location, if the melt continue to flow after mending, and the quality of mending based on the material's temperature, pressure and level of mixing (particles of different colors blended together).

[0152] At 6K, the method checks for the time. Stop time is a user defined variable. If stop time has not been reach, returns the module to 6G for simulation at a different time step. If it is stop time, the simulation is completed at 6L. User may output result of the simulation.

[0153] While the invention has been illustrated and described in detail in the drawings and foregoing description, the same is to be considered as illustrative and not restrictive in character, it being understood that only the preferred embodiment has been shown and described and that all changes and modifications that come within the spirit of the invention are desired to be protected.