Simulation Of Micro Injection Moulding With Emphasi s On The Formulation Of Feedstock Viscosity : Use Of Non-equilibrium Molecular Dynamics For The Determination Of Viscosity Of Multi-body Fluid

The need for prediction of shear viscosity of fluid in particle charged Micro Injection Molding at mesoscale, by modelling a whole system p article-polymer with inter-dependencies, permits to establish a more realistic feedstock viscosity f ormulation. The applicability of non-equilibrium molecular dynamics (NEMD) is investigated for the d termination of shear viscosity of melts composed of particles/polymers in microcavities. NE MD is used to simulate planar Poiseuille flow of metallic particle-polymer melt. Simulations are car ried out using molecular dynamics simulation package ESPResSo. The variation of viscosity face t o t mperature is in agreement with theoretical results. Simulations are compared to experiments. T he equivalent viscosity formulation is tuned according to NEMD simulation results, and implement ed in a MIM solver built up by the authors. MIM simulations are compared to previously implemen ted simulations using another equivalent viscosity formulation based on experiments for the case of mono-injection moulding. Résumé. Le papier propose une approche de la viscosité à l’échelle mé soscopique de cisaillement, pour des mélanges fortement chargés en poudres util i és en moulage par injection de poudres métalliques. Cette modélisation pour la prédiction plus réaliste de la viscosité est basée sur l’étude des interactions d’un système liants thermoplastiqu es/particules en utilisant une méthode de dynamique particulaire (DPD) en particulier la méth ode NEMD (non-equilibrium molecular dynamics). Les développements numériques ont été ré alisés en utilisant le logiciel ESPResSo. Des simulations de remplissage de micro-cavités ont été r alisées en prenant en compte la viscosité à l’échelle mésoscopique et un modèle équivalent a ét é d veloppé par les auteurs pour prendre en compte la taille des particules ainsi que le taux d e charge en poudres. Les résultats numériques obtenus sont prometteurs et encourageants pour l’av enir.


Introduction
The knowledge of shear viscosity in Metal Injection Moulding (MIM) is a key factor to carry out realistic numerical simulations [6,7].So far, the viscosity of fluids is determined by the mean of experimentations associated to phenomenological laws that involve several parameters such as polymer grades, powder contents, metallic particle size and temperature of the melt.Statistical mechanics provides a well-established link between microscopic equilibrium states and thermodynamics.If one considers systems out of equilibrium, the link between microscopic dynamical properties and non-equilibrium macroscopic states is more difficult to establish [1,2].For systems lying near equilibrium, linear response theory provides a way to derive linear macroscopic laws and the microscopic expressions for the transport properties that enter the constitutive relations in these laws.If the system is displaced far from equilibrium, no fully general theory exists to treat such systems.By restricting consideration to a class of non-equilibrium states, arising from perturbations (linear or non-linear) of equilibrium state, methods can be developed to treat non-equilibrium states.Furthermore, non-equilibrium molecular dynamics (NEMD) simulation methods can be devised to provide approximations for transport properties of these systems.In non-equilibrium molecular dynamics simulations of steady states, the thermostat is of great importance.The thermostat used here is the method of Dissipative Particle Dynamics (DPD).Many outlines of DPD discuss two aspects, first, soft particles which should represent a cluster of atoms, second, a momentum-conserving stochastic thermostat added in order to model the internal degrees of freedom, which results in dissipation.As it is stipulated in [3], it is legitimate to use the DPD thermostat also for simulations with hard particles.A standard Molecular Dynamics (MD) system with an added DPD thermostat is run, thereby being able to afford a substantially larger time step compared to pure MD, and nevertheless correctly reproducing hydrodynamic behaviour.

Feedstock viscosity formulation
Since the viscosity relationship established by Einstein, many studies have attempted to establish phenomenological equations to model the viscosity of a concentrated medium.The first of them have emerged in the mid 20 th century.A new model, named equivalent viscosity model, is proposed to take into account 4 parameters: solid volume fraction φ , shear rate γ& , temperature T and particles size d 10 , d 50 and d 90 .This model is a power law viscous terms, an advection energy term, a term that is related to solid loading and a term that account to the powder size in an average way.The viscosity value is expressed as: where η is the viscosity, η 0 the apparent viscosity, φ the solid volume fraction, φ m the maximum solid volume fraction, d 10 particle size at 10%, d 50 particle size at 50%, d 90 particle size at 90%.Parameters A and B are determined to fit best the experimental data.The feedstock is modelled as a system composed of FeNi metallic particles and polymeric binder.

NEMD Algorithm for Shear
Green and Kubo [9] showed that the phenomenological coefficients describing many transport processes and time-dependent phenomena in general could be written as integrals over a certain type of function called a time-correlation function.The Green-Kubo formulas are the formal expressions for hydrodynamic field variables and some of the thermodynamic properties in terms of the microscopic variables of an N-particle system.The identification of microscopic expressions for macroscopic variables is made by a process of comparison of the balance and conservation equations of hydrodynamics with the microscopic equations of change for conserved densities.The importance of these formulas is three-fold: they provide an obvious method for calculating transport coefficients using computer simulation, a convenient starting point for constructing analytic theories for non-equilibrium processes, and an essential information for designing nonequilibrium molecular dynamics algorithm.The Green-Kubo formula for the shear viscosity is given by where P αβ is an off-diagonal (α ≠ β ) is expressing the viscous pressure tensor.k is the Planck constant.V and T are the volume and temperature of the polyatomic fluid, respectively.Angle brackets denote an equilibrium ensemble average.
where v iα and v iβ are the α and β components of the velocity of particle i, r ijα is the α component of the distance between the particles i and j, and F ijβ is the β component of the force of their interaction.
The general principle of the Non-Equilibrium Molecular Dynamics (NEMD) method [1] is to introduce a possibly fictitious external field X into the system motion equations which derives the corresponding thermodynamic flux J .The first requirement for this applied field is that it should be consistent with the periodic boundary conditions to ensure that the simulation box remains homogeneous.The second requirement is that the transport coefficient γ of interest can be calculated from the constitutive relation: Historically two approaches were proposed to simulate fluids under shear: the boundary driven algorithm and the Sllod algorithm.In the boundary driven method, a modification of the boundary conditions suffices to simulate systems far from equilibrium.For example, Lees-Edwards boundary conditions [8] can be used to model shear flow.The shortcoming of this approach is the lack of an explicit external field as driving force in the equations of motion.Therefore, it is difficult to link the results with the Green-Kubo relations or the study transitional flows.
The Sllod algorithm [1], a standard method using homogeneous Lees-Edwards sliding brick boundary conditions, still uses Lee-Edwards boundary conditions, but also incorporates the flow field into the equations of motion.It is the most efficient technique for calculating the shear viscosity.This method has been chosen to study polymeric chains and metallic particles under shear.This algorithm sets up a steady state planar Couette flow with the two plates moving in opposite x directions located at y = ±∞ so that the streaming velocity has a non-zero component in the x direction.
Following Evans and Morriss [1] the equations of translational motion for the center of mass in a molecular fluid are given by: where u = u x ,0,0

(
) with u x = γy is the velocity field corresponding to planar Couette flow.r i , m i and p i are respectively the position vector, mass and momentum of particle i. F i is the force applied to particle i.These equations of motion are combined with the Lees-Edwards "sliding brick" boundary conditions [8].In the absence of the thermostat and the isobaric constraint, the terms in equation ( 4) involving the strain field, γ , cancel to yield Newton's equations of motion relating r i and F i .This implies that the Sllod algorithm truly generates boundary driven planar Couette flow, leading to the conclusion that it is correct to arbitrary order in the strain rate [1].In order to obtain a good signal-to-noise ratio, with NEMD it is necessary to use strain rates γ which are high enough to cause the shear viscosity to be strain rate dependent.In order to compute the shear viscosity of a Newtonian fluid using the Sllod algorithm, after the simulation reaches the steady state at a given strain rate γ one computes and averages the pressure tensor defined in equation ( 2).The strain rate dependent shear viscosity is then obtained from Newton's law of viscosity where P xy and P yx are the averaged xy and yx components of P .From kinetic and mode coupling theories, to leading order the strain rate dependence of the shear viscosity is linear in γ 1 2 .Hence, to apply the Sllod algorithm to a Newtonian fluid, one performs several simulations at differing strain rates γ and fits the resulting strain dependent viscosities to the equation: The zero strain rate extrapolation of η, η 0 , corresponds to the Newtonian viscosity.

Shear Viscosity Simulation Setup
The hydrodynamics boundary conditions result from interactions between the fluid particles and the walls.Depending on the microscopic structure of the wall/fluid interface, these can be quite complex.In this approach, boundaries are replaced by hard planar surfaces, and the unknown atomistic forces by an effective coarse-grained friction force between the fluid particles and the walls.This leads to a dissipation of the kinetic energy and therefore to a deceleration velocity of the fluid close to the boundaries.The resulting slip length depends on the strength of the friction force.The effect of boundary conditions is modelled using partial-slip boundary interactions added to Lennard-Jones boundary interaction to mimic the wall repulsive effect.
The system is composed of N particles and M polymers composed of MB beads each.A squared-section channel of size L x *L y *L z is modelled.Four thin layers perpendicular to the y axis, with distance L y /4, were chosen for driving in the +x and -x direction, respectively.
Interparticle spherically symmetric potentials are continuous in order to set up a standard Molecular Dynamics procedure using the Velocity-Verlet algorithm [4], and short-ranged, in order to keep the number of force calculations at minimum.A convenient choice for this is a truncated and shifted Lennard-Jones 12-6 potential.One will henceforth use the Lennard-Jones units where LJ parameters are chosen according to material characteristics.Beads of three different species (Polypropylene PP, paraffin wax PW and stearic acid SA) are considered, which interact through bonded interactions.From the polymer simulations it is known that it is computationally efficient to link the dimmers via harmonic Finite Extensible Nonlinear Elastic springs.After experimental verification in laboratory, Coulomb interactions are neglected faced to other interactions.
Figure 1(a) shows shear viscosity simulation results obtained using ESPResSo [5] for a Poiseuille flow in a square section micro-channel for a melt dense at 80%, at several temperatures k B T (noted T in Figure 1).

Simulation of metal injection molding injection stage
Numerical simulations of the injection process were carried out for 2D mold cavities with the implementation of the equivalent viscosity formulation.In case A, a Finite Element Method is employed to solve Navier-Stokes equations and Streamline-Upwind/Petrov-Galerkin method for the solving of the transport equation [6]. Figure 2 relates the advance of the filling front for both constant ( µ = 10 3 Pa.s ) and equivalent viscosity formulation defined by equation (1).Global filling ratio progress is highlighted in figure 3. No-slip boundary conditions are imposed, velocity is imposed on inlet and a zero pressure imposed on outlet.In case B, the solving procedure involves a combination of Finite Differences/Finite Element Methods for solving the coupled Navier-Stokes and transport equations [7].Comparisons of simulations implemented with both constant ( µ = 10 3 Pa.s ) and variable viscosity formulation (equation ( 1)) are highlighted.Figure 4 stands for the advance of the filling front for both constant and variable viscosity formulation.Global filling ratio progress is highlighted figure 5. No-slip boundary conditions are imposed, velocity is imposed on inlet and traction condition is imposed on outlet.For both cases, one will note that less diffusion of the filling front and more complete filling of the mould cavity are achieved using the variable viscosity formulation.

Conclusion
NEMD is revealed a proper tool to describe a particle charged fluid under shear strain rates to obtain the shear viscosity.Shear viscosity results are in agreement with the theory concerning temperature variations.Numerical simulations of the injection molding process show very interesting results by the use of the proposed equivalent viscosity formulation.A further step in the implementation of NEMD simulations using ESPResSo will consist in simulating feedstocks with different powder compositions and binder/particle ratio in order to tune the equivalent viscosity formulation.

Figure 1 (
b) exposes results for two different feedstocks with different powder charge loading and binder proportions.As expected, a temperature raise implies increase of shear viscosity.High proportion in powder charge implies decrease of shear viscosity.Powder charge loading has a more obvious impact on shear viscosity than DPD temperature.

FIG. 1 .
FIG. 1.(a) Shear viscosity of feedstock composed of 60 vol% metallic particles and 40 vol% binder at several DPDtemperatures, (b) Shear viscosity of feedstocks at DPD temperature 1 for several feedstocks composition.

FIG. 2 .FIG. 3 .
FIG. 2.Filling front advance at for significant time steps for constant (a) and equivalent viscosity (b) formulation in the case A.

FIG. 4 . 6 FIG. 5 .
FIG. 4. Advance of the filling front at for significant time steps for constant (a) and equivalent viscosity (b) formulation in the case B.