A simple and systematic scheme implemented in explicit FEM solver for surface tension effects in powder injection moulding process

In order to simulate more accurately the powder injection moulding process, the explicit finite element method solver is extended with the surface tension effect. The evaluation of surface tension takes the notion of pressure boundary method, while a simple and systematic scheme is proposed to fit the finite element method solver for the Laplacian operator. Because of the difference in dimension for filling function and velocity function, the integration of filling function in second derivative is not suitable to be transformed into the boundary integration and the integration of function in lower order derivative. To evaluate conveniently the curvature of filling front, hence the force of surface tension, a simple and systematic scheme is suggested and implemented into the finite element method solver. This specific scheme includes only the vectorial operations in low cost, and is completely systematic without piecemeal operations. Fitness of the proposed method is proved by the numerical examples of filling flow in a small-scaled channel. It shows the considerable effect of surface tension for the problems in micro-scale of sub-millimeter sizes, in which the boundary conditions at front surface are not negligible in powder injection moulding process. The surface tension effect becomes the dominating role for governing the trace and shape of filling front, which can no longer be neglected.


Introduction
Injection moulding process is an effective manufacturing method for mass production of the components with very small size and complicated shape [1]. Due to size decrease, filling flow in moulds under a certain scale exhibits the rheological phenomena different from the conventional ones [2]. Some factors as surface tension, wall slip, viscous dissipation and convective heat transfer will play the important roles in small-scaled injection, even though they can be ignored in the analysis of large-scaled ones. In order to evaluate the surface tension effect of injection flow in the small-scaled die cavity mould, a special scheme is proposed and implemented in explicit Finite Element Method solver [3] to fit the Pressure Boundary Method (PBM). Then the surface tension effect in injection moulding process is evaluated.
In injection moulding process, thermoplastic polymer is injected into the die cavity mould then air inside is squeezed out. In this process, there exists a moving interface between the polymer and the air. Because of the large difference in material properties, molecules at the interface are significantly more attracted by the molecules of polymer than that of the air. Then the molecules at the interface subject to the resultant force directed inward of the polymer.
The surface tension effect in micro-injection has been previously studied by W. Cao et al. [4] and Kim et al. [5]. W. Cao et al. developed a model to take into account the dynamic contact angle between fluid and micro die cavity to simulate the surface tension effect in the filling of micro-channel. Kim  In the paper, to evaluate the curvature of filling front, hence the force of surface tension, a simple modelling scheme is proposed and implemented into explicit finite element software. The filling ratios in mould die cavity at the same instant, with or without consideration of surface tension effect, will be compared. The comparisons are made for both the sizes in conventional injection moulding and the sizes for microinjection process. For the investigation of powder injection moulding process, some simulations are made for the example of filling flow in sub-millimeter sizes of high loaded polymer in a straight channel of features heaving.

Methods for capturing interface and curvature
The simulation of injection moulding at micro-scale is dependent on accuracy and stability of the surface tension evaluation [6]. The feedstock of filling flow may consist of high loaded polymer with metallic powders. The effect of surface tension force in filling flow is shown in Fig. 1. Front capturing is prerequisite for implementation of the proposed scheme.
Some numerical methods for front capturing had been developed over the past few decades. The well-known methods are volume-of-fluid [7], level-set [8], phase field [9], and smooth volume of fluid [10] methods. Among them, the coupled volume-of-fluid (VOF) and level set (LS) methods (CLSVOF method) have been extensively employed. The CLSVOF method was introduced by Bourlioux [11], and then put forward by Sussman and Puckett [12], Son and Hur [13] and Menard et al. [14]. By virtue of the CLSVOF method, curvature of the front surface can be computed accurately. The interface is reconstructed via a piecewise linear interface construction (PLIC) scheme from the VOF function. While the normal vector of interface is computed from the LS function n ! ¼ ∇φ= ∇φ j j (where φ is a function of the signed distance to interface). The principle of this coupling is to use VOF function for mass conservation, while the interface and curvature is captured smoothly by the LS function [15]. Curvature of the interface is evaluated by: where n ! is the interface normal vector, φ is an approximation of the signed distance function to the interface. To solve the differential (Eq. 1) for curvature, the approaches such as the Finite Volume Method (FVM), Finite Difference Method (FDM) and Finite Element Method (FEM) are widely used. It is discretized by the finite volume method on the non-staggered meshes in [16]. Sun and Tao [17] proposed an iterative geometric operation to calculate the level set function near interfaces, which can be applied to compute the accurate curvature and smooth the discontinuous physical quantities near interfaces. Based on FVM, it uses the high order essential non-oscillatory (ENO) upwind difference [18] for the convection term and the central difference for the viscous and curvature terms. In [19], the advection equation of level set function is solved by a second-order finite volume method. The advection of volume fraction is performed by a scheme based on bounded compressive normalized variable diagram (NVD). The interface is reconstructed by both the level set and the volume fraction information. In [12], the CLSVOF method is used but it does not smooth the curvature at all. The curvature is obtained via finite differences of the level set function, which is derived at the previous time step by both the level set function and volume-of-fluid function.
Compared to finite difference method, finite element method is easier to be used with a great variety of element types, and is thus capable of handling complex geometries and boundary conditions. Hysing [20] introduced a new level set method for simulation of immiscible fluid flows, which essentially combines a non-conforming finite element flow solver with a conforming level set interface tracking method. It uses a more accurate approach to reconstruct the curvature via L2-projection and patch recovery techniques. In [21], surfaces are implicitly represented by the level set method. It studies the possibility to use the shape functions of higher order approximation with extended finite element method (X-FEM), when geometries include curved discontinuities. The finite element approximation is enriched by additional functions with the notion of partition of unity to track material interfaces. Yang et al. [22] presents an adaptive CLSVOF method for interfacial flow simulations on unstructured triangular grids. At each time step, it evolves both the level set function and the volume fraction. Evolution of the level set function is governed by advection equation. A discontinuous Galerkin finite element method is used for its solution. The advection of volume fraction is performed by a Lagrangian-Eulerian method. Methods for solving the Laplacian operator As a key issue to implement the surface tension in finite element methods, the curvature of filling front should be evaluated as illustrated in (Eq. 1). It represents a second derivative of the function, in the form of a Laplacian operator. As the filling function is a scalar one while the velocity function is a vector one, the integration of filling function in second derivative is not suitable be transformed into the boundary integration and the integration of function in lower order derivative. It is difficult to use ordinary integral transform to reduce the derivative order of integrated function. Some methods were previously proposed to settle this problem. In traditional level set methods [23], curvature is computed from the spatial derivatives of a scalar function at any instant in time, as κ (where φ is a function of the signed distance to interface). On a Cartesian grid, this divergence is generally calculated at node points by using a nine-point stencil with center differences for all the partial derivatives. Ramšak and Škerget [24] had developed an efficient 3D multi-domain boundary element method (BEM) for solving problems governed by the Laplace equation, where the domain problem is transformed into a boundary problem. It uses a multi-domain approach, also known as the subdomain technique, very similar to FEM. Trontin et al. [25] used markers with a particle/ level-set method to redefine the curvature between two immiscible phases. The cubic spline interpolation gives each marker an estimation of the local curvature of the interface. It can be evaluated exactly at the interface without any interpolation from the Eulerian grid. Rangogni [26] deals with the generalized Laplace equation by coupling the boundary element method with the perturbation method. The method requires an internal network but the unknowns are only on the boundary. Based on the VOF method, Meier et al. [27] devises a new method, which uses empirical formulas obtained from databases that has been generated and stored in a data bank to determine the interface curvatures. Raessi et al. [28] presents a new method for calculating interface normal vectors and curvatures, where the interface unit normal is advected along with whatever function represents the interface, and curvatures are calculated directly from these advected normal.

Original work of the present paper
Due to the lack of appropriated FEM method for curvature calculation, the present investigation proposes a simpler and easier way for implementation of the surface tension effect in finite element method. The work provides an effective method for the analysis of capillaritydominant free surface flows in small-scaled injection moulding, which involves the computation of Laplacian operation for evaluating the curvature of filling front. It can be used in the problems that require the integration of a scalar function (filling function) in second derivative, but the effect is associated with the vector functions. It affords the practical way for evaluation of the surface tension effects in viscous filling flow. The incorporation of surface tension effect is realized under the frame of Taylor-Galerkin method, which is used for solution of the advection equation to track the filling front. The notion of CLSVOF method is used. The force of surface tension is incorporated into the Navier-Stokes equations via an individual source term. At each time step, both the level set function and the volume fraction are evolved to update the filling front. The level set function in advection equation is evolved by the Taylor-Galerkin method of second-order characteristics. The interface is reconstructed based on the information of level set and volume fraction. It provides the necessary accuracy for evaluation of the surface curvature and the surface normal on the filling front.
The present work takes the mathematic notion of surface tension force as explained by Tong and Wang [29] with the easy method proposed and the implementation is realized in our finite element software [3]. The work of Tong and Wang [29] is based on the finite difference method with the easy way to evaluate the second derivative of level set function. The FEM explicit solver employs the simple elements of low order interpolation functions, so no second derivative of the function can be obtained by a direct way. The numerical implementation is clearly distinguished. The explicit finite element solver, developed by the research team, appears to be an efficient tool for simulation of the viscous filling flow [3]. However, the elements of low order interpolation functions do not provide the direct way to evaluate the second derivative of represented function. To avoid the complication in evaluation of the Laplacian operation on filling function, a specific method is proposed to evaluate curvature of the filling front by the simple and systematically procedures. Then the force of surface tension can be introduced directly in solution of the Navier-Stokes equations. By comparison between the magnitudes of surface tension force with the viscous force for a filling example in small-scaled channels, the effects of surface tension in injection moulding process are studied. It shows the importance of surface tension in sufficiently small-scaled injection moulding of sub-millimeter sizes, though this effect does not represent the significance in ordinary injection moulding process.

General definition
As usually chosen, Eulerian description is adopted for simulation of the mould filling problems, which avoids the complicated and expensive remeshing procedures with Lagrangian description. The general definition of mould filling problems with surface tension is expressed as followings: Let t ∈ [0, t 1 ] be an instant in the injection course, in which t 1 is the last moment of filling process to reach the fully filled state. The sum of position X in the whole mould die cavity is defined as set Ω. The set Ω consists of two different portions at each instant, the portion Ω F filled by thermoplastic polymer and the remained space Ω V taken by the air. Based on the VOF definition, a field function F(x,t) is defined to represent filling state of the injection moulding at different instants. This field function takes value 1 to indicate the portion filled by thermoplastic polymer and value 0 for the remained void portion, which contains in fact the atmosphere. The physical and geometrical definition for this modelling is shown in Fig. 2. In which Γ I indicates inlet of the die cavity mould, Γ O represents the outlet through it the air originally in the die cavity mould can be squeezed out during the injection. Γ S stands for intersection of the subsets Ω F and Ω V , which is in fact filling front of the injection flow. The resultant of surface tension is produced on the surface of filling front and exerted towards center of the curvature inside the filling material. P is the injected pressure imposed on inlet boundary. P o is the pressure on outlet boundary, which is set to be 0 in the present work to represent the environment pressure. According to the continuum surface force (CSF) model [30], the surface tension force F st is handled to be a body force F b in a thin layer of the front. The details will be discussed later in the following contents.

Governing equation for filling process and evolution of filling front
To keep singleness of the solution strategy and simplicity of the software structure, under the frame of Eulerian description, the governing equations for filling flow are chosen the same for filled and void portion of the injection model, except that the physical parameters are chosen differently for these two domains. In fact, the interest of simulation is the filling flow of high loaded polymer in die cavity.
The front position of filled domain is represented by the predefined filling function F(x,t). This variable is governed by an advection equation, driven by the velocity field in die mould cavity.
The boundary condition is F=1 on inlet of the die mould. Its initial condition is F=0 everywhere in the die mould except for the inlet surface. For the convenience in numerical methods, the solution of above advection equation is based on the continuous functions, and discretized at each time step. The continued function is in fact the LS function, named φ by other authors [31]. It is called the filling function F, for clearness and integrity of the presentation. It changes continuously across the filling front, with distribution from higher value in filled portion to lower value in void portion. The solution of advection equation (Eq. 2) is to determine the function F(x,t), with solution procedure same as the previous work in research team [32]. Then the filling front locates at the position X F where F takes a prescribed value between the larger value in filled portion and smaller value in void portion.
where F P is a prescribed value to identify the position of filling front according to the distribution of filling function F. This value can be adjusted for the purpose of mass conservation in filling process.

Momentum conservation of filling flow
The continuum surface force (CSF) model [30] has been widely used to introduce the surface tension. In the CSF model, surface tension effect is treated as a body force F b in a thin layer of the elements that locate on the filling front. It can be added directly into the Navier-Stokes equations as a source term. Considering surface tension force, the equations of momentum conservation for two different portions can be expressed as: where ρ p is the density of feedstock in the filled die mould cavity, ρ a is the air density in the unfilled die mould cavity. P represents the hydraulic pressure field, σ 0 p and σ 0 a are the deviatoric Cauchy stress tensors in filled and void portion, that is, in the feedstock and air inside the die mould cavity. g is the gravity vector. F b represents the surface tension force in the filling front. It should be mentioned that the material properties for air portion are chosen different from their true values, for the purpose to keep numerical stability. It should be mentioned that the material properties for air portion are chosen different from their true values, for the purpose to keep numerical stability. In the vectorial algorithm with finite element method proposed by Cheng [32], the solution scheme in air portion and filled portion is the same to keep singleness of the numerical strategy. It is necessary and effective for software realization except for a small problem. It may result in instability of the numerical solution if the material properties in void portion (occupied by air in fact) and filled portion are too much different. For the sake of numerical stability, the material in the void portion is modified some different from real properties of the air. In fact, simulation result in the void portion is not of our interest. Thus, such a numerical treatment is acceptable.

Surface tension force
With the CSF model, surface tension effect is treated as a body force F b and added into the Navier-Stokes equations as an external force [30,33]. It is distributed within a transition region of finite thickness at the interface, given by: where σ is the coefficient of surface tension, κ is the curvature of filling front surface, n ! is normal vector of the surface, indicated inward the injected material. δ(x) is a delta function concentrated at the interface. From the definition of filling function F(x,t), the gradient of filling function ∇F(x, t) represents a vector in inner normal direction of the filling front. Its value in tangential direction of the interface is zero. So the following expression can be used: The above expression serves to determine unit normal of the filling front surface, no matter how much is the value of filling function and the magnitude of its gradient.
The body force in (Eq. 6) can also be expressed into: This body force will be included into the momentum equation as a source term, as F b in (Eq. 4). This continuous treatment of the discontinuous change at the interface makes easier the implementation of surface tension effect. In problems with complex topological changes, the CSF model is superior to the conventional method in robustness and versatility.

Implementation of surface tension in finite element method
From (Eq. 8), it is evident that accuracy of the curvature estimation is the key issue in implementation of surface tension.

Surface curvature computation
As in (Eq. 1), normal and curvature of the filling front can be calculated by the distance functions in LS method. In mathematics, it depends on the shape of filling front surface, no matter value of the LS function and the magnitude of its gradient. In the present study, the LS function is replaced by the filling function F(x,t). The expression in (Eq. 1) then becomes: The filling function F(x,t) is determined by the solution of advection equation (Eq. 2). Based on the filling function F at actual instant, curvature of the filling front κ(x) can be determined by an explicit way.
For implementation of the surface tension in finite element method, it is necessary to realize the second derivative of filling function, incorporating with the interpolation functions of finite elements. The previous research resulted in the development of an explicit algorithm, in which the simple elements of low order interpolation functions are used [32]. The integration of filling function in second derivative is not suitable be transformed into the boundary integration and the integration of function in lower order derivative, because of the difference in dimension for filling function and velocity function. To evaluate conveniently the curvature of filling front, a Fig. 3 Finite element mesh of straight channel for simulation of injection filling process in 2D case simple and systematic scheme is suggested and implemented into the explicit FEM solver. This specific scheme includes only the vectorial operations in low cost, and is completely systematic without piecemeal operations. In solution of the Navier-Stokes equations, the source term including field κ (x) is determined by the following procedure: 1) First, obtain the gradient field at the nodes where A represents the assembling operation of finite element methods, [N] is the matrix of interpolation functions for filling function.
[G] is a gradient operator consists of the derivatives of interpolation functions for filling function F. {∇F} stands for the discretized field of ∇F, represented at the nodes. {F e } includes the node values of F for a prescribed element.
2) In the second steps, the Laplace of F is determined as follow The values determined at Gauss point can be obtained by: If the discretized values at nodes ΔF are required for next procedure of the solution, the discretized field of Laplace F is determined by where [D] is a divergence operator, built based on the interpolation function for filling function.

Surface tension force in computation
In finite element method, the surface tension is implemented in a layer of the elements, located at the positions of filling front surface. It is processed as body force acting in a thin layer at the front. Term F b in (Eq. 6) can be directly added into the Navier-Stokes equations, as shown in (Eq. 4). Based on the explicit algorithms, the evolution of velocity field is effected by two fractional steps [3]. Effects of viscous diffusion, force loading and surface tension are considered in the first step. The incompressibility is kept in the second step, with the solution of pressure field. By Galerkin method [34], Navier-Stokes equations (Eq. 4), including the term of surface tension, can be discretized in the following form: where V * n is the velocity fields of temporary middle values, M is the mass matrix lumped into diagonal form, F σ is the viscous diffusion term, F ext is the load vector with gravity contribution and boundary conditions. F st is the term of surface tension force.
where A represents the assembling operator in finite element method, B is the matrix of derivatives of the interpolation functions [N] in an element, defined to calculate the symmetric part of velocity's gradient ε ¼ BV e , V e represents the nodal values of velocity vectors in an element, σ ′ stands for the deviator of Cauchy stress tensor, F ext stands for the surface loads converted to the nodes on boundary.
The second fractional step is to satisfy incompressibility of the filling flow, with the solution of pressure field. The relevant literatures can be referred for its implementation [3].

Numerical investigation and discussion
In order to study the effect of surface tension in injection moulding process, the evaluation of surface tension force in filling viscous flow is realized by simulation. The development is made on the explicit algorithm with fully vectorial operations [3], previously developed in research team. Based on the finite element software developed on Matlab platform, this section focuses on realization and validation of the new developed function for surface tension. elements is shown in Fig. 3. The width of mould inlet W = 10 mm and the mould length L = 5 × W = 50 mm. 316 L stainless steel feedstock density, corresponding to 316 L stainless steel powders and thermoplastic binder based on paraffin wax, is set to be 711 kg.m −3 . Feedstock shear viscosity is 20 Pa.s. A constant pressure P = 10 MPa along the axial direction is imposed on mould inlet. Surface tension coefficient is supposed to be a constant σ=0.5. These material properties and boundary conditions are assigned just for the purpose of software validation. They do not mean the real values obtained from experiments. As the present development and analysis is based on the explicit algorithm of research team [3], some technical issues are quoted directly. The present work focuses on investigation of the surface tension effects in filling processes of the viscous flow. The comparison of filling states with or without the effect of surface tension at the same filling instant is shown in Fig. 4. The component of surface tension force in horizontal direction x and in vertical direction y is shown in Fig. 5.
For analysing the impacts of surface tension on mould filling processes of different scales, the examples of different sizes are realized in simulation. In order to understand the variation of surface tension effect when the geometrical scale of mould filling changes, the examples are made under the same conditions, except for their sizes. The meshes for different models of different sizes remain proportionally the same, except that the size of elements changes accordantly. Then the influence of element distribution can be removed for the examples of different sizes. Figure 6 (including 6.1~6.6) represents variation of the simulation results for filling state at different geometrical scales, with disregard and consideration of the surface tension term. The comparison of two results is labelled by (a) and (b) for the simulations of each size, respectively. The notation W and L are shown in Fig. 3. W is the width of mould die cavity inlet, and L = 5 W is mould die cavity inlet length.
However, results in 6.3~6.6 seems to be a little over a real filling patterns. It's caused by the unrealistic material properties and boundary conditions, which are not the precisely one obtained from experiments, they are assigned just for the purpose of software validation. So the evolution of filling front represents only the referential shape. It is certain that it will make no sense to the conclusion of surface tension effects. The really reliable parameters of injection moulding flow should be achieved by the experimental ways in the further works.
Comparison of the simulation results for different sizes, shown in Fig. 6, provides variation of the surface tension effect according to the change in geometrical scale. All simulation examples have been realized with the same material properties and the shear viscosity behaviour for the flow of feedstock. The difference in curvature of the filling front surface is observed. The smaller the example is, the more obvious for the difference in surface curvatures. It indicates that the surface tension force impacts on the shape of filling front in injection moulding process. It increases the curvature value of filling front. These results explain well the physics nature of surface tension, which tells that the surface tension induces contraction of the liquid surface. Meanwhile, it should be noticed that the evolutions of filling state with and without the surface tension term are nearly the same while the size scale is larger than 1 mm. The surface tension effect becomes more important and cannot be neglected in micro scaled injection moulding when channel width is within 1 mm. The filling ratio in mould die cavity at the same instant, with or without consideration of the surface tension effect, is shown in Table 1. With F filling ratio in mould die cavity without surface tension and F 1 filling ratio in mould die cavity with surface tension The data in Table 1 indicates that the surface tension contribution to micro cavity (sub-millimeter sizes) is considerably larger than that to macro cavity (above millimeter sizes). In addition, the surface tension takes negative effect on the filling process. It may hinder the filling evolution. The value of surface tension force increases obviously when the mould die cavity size decreased. This phenomenon can be also found in Table 2 and Fig. 7.
For the flow of small Reynolds number, filling front surface is governed by advection equation, advected by velocity field. The solution of velocity field is affected by viscous force and surface tension force. In order to evaluate and compare their effects on advance of the filling front, 10 nodes on filling front surface were chosen to compare the values of surface tension and viscous force. The average value of surface tension force and viscous force at 10 nodes was used to evaluate the effect of surface tension on evolution of the filling front surface, as shown in Table 3.
From Tables 2, 3 and Fig. 7, it can be found that neither viscous force and surface tension force is the only driven mode during filling period. The surface tension effect increases obviously when the mould die cavity size decreases. Especially, for the sub-millimeter sizes case in micro cavity, the average value of surface tension force represents more and more close to or even larger than the viscous force. The Weber Number also reduces sharply close to or even less than 1.0. It indicates that surface tension force will play the dominating role for governing the trace and shape of filling front. This conclusion agrees well with the work of Chien-Te Li et al. [35], who investigated the effects of gravitation and surface tension on the evolution of flow front in injection of the center-gated disks.
To validate the conclusions above, other three models with more complex cavities are specially designed for the simulation, as shown in Fig. 8. All of them have variation of the cross sections along the flow paths. In order to learn the surface tension effect on evolution of the filling state in each of the models, two similar cases are used to make the comparison work. All conditions are set to be the same for these two cases except the consideration or not the surface tension effect. The same pressure is imposed on the inlet boundary. The surface tension effects can be observed by the difference between two filled ratios versus time, as shown in Fig. 8.
From Fig. 8, the variational trend of surface tension effect according to the mould die cavity sizes is obvious. All of the three mould die cavities have the variation of cross section from a bigger one to the smaller branches. It makes accordingly the difference of filled ratios at different instant for the cases with or without surface tension effect. The difference in filling ratios of these two cases is growing from a very low value to the peak value until one of them getting to the fully filling state, then the other case continues to fill until difference of the filling ratios become zero.

Conclusion
The work describes a novel and simple methodology to implement the surface tension in FEM solver. It provides an effective way to evaluate the complex effect of surface tension in injection moulding of small powder injection components. In micro injection moulding or polymer processing at the sub-mm scale, the tension at front surface represents the important effect on filling of the feedstock. The effect of surface tension has been implemented in an explicit finite element solver self-developed in laboratory. This software is dedicated to the filling flow of high loaded polymer for powder injection moulding process. A systematic scheme has been proposed and realized to fit the requirement of finite element methods. The Laplacian operation of filling function has been realized by a twostep vectorial process to evaluate curvature of the filling front. The results are introduced directly into the source term of surface tension force in the solution of Navier-Stokes equations. According to the continuum surface force model, the force of surface tension is distributed in a thin layer of the elements, which locates on the filling front of injection flow. term increases when the mould die cavity size decreases. The paper shows validity of the model at mm-scale and sub-mm scale. These numerical results are in close agreement with the data reported by Li and al. in [35]. The surface tension can be neglected in conventional macro injection moulding processes above millimeter sizes, but it represents significant effect in filling of the micro die cavities with sub-millimeter sizes.