# Simulation of micro injection moulding with emphasis on the formulation of feedstock viscosity: use of non-equilibrium molecular dynamics for the determination of viscosity of multi-body fluid

Guillaume Larsen, Zhi Qiang Cheng, Thierry Barrière, Bao Sheng Liu, Jean-Claude Gelin

<figure><img src="image_1.png" /></figure>

## To cite this version:

Guillaume Larsen, Zhi Qiang Cheng, Thierry Barrière, Bao Sheng Liu, Jean-Claude Gelin. Simulation of micro injection moulding with emphasis on the formulation of feedstock viscosity: use of non-equilibrium molecular dynamics for the determination of viscosity of multi-body fluid. International ESAFORM Conference on Material Forming, Apr 2011, Belfast, Ireland. 10.1063/1.3589599. hal-02300610

HAL Id: hal-02300610

https://hal.science/hal-02300610

Submitted on 22 Feb 2023

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Simulation Of Micro Injection Moulding With Emphasis On The Formulation Of Feedstock Viscosity: Use Of Non-equilibrium Molecular Dynamics For The Determination Of Viscosity Of Multi-body Fluid

Guillaume Larsen<sup>a, b</sup>, Zhi Qiang Cheng<sup>b</sup>, Thierry Barriere<sup>a</sup>, Bao Sheng Liu<sup>b</sup> and Jean-Claude Gelin<sup>a</sup>

**Abstract.** The need for prediction of shear viscosity of fluid in particle charged-Micro Injection Molding at mesoscale, by modelling a whole system particle-polymer with inter-dependencies, permits to establish a more realistic feedstock viscosity formulation. The applicability of non-equilibrium molecular dynamics (NEMD) is investigated for the determination of shear viscosity of melts composed of particles/polymers in microcavities. NEMD is used to simulate planar Poiseuille flow of metallic particle-polymer melt. Simulations are carried out using molecular dynamics simulation package ESPResSo. The variation of viscosity face to temperature is in agreement with theoretical results. Simulations are compared to experiments. The equivalent viscosity formulation is tuned according to NEMD simulation results, and implemented in a MIM solver built up by the authors. MIM simulations are compared to previously implemented simulations using another equivalent viscosity formulation based on experiments for the case of mono-injection moulding.

**Keywords:** Injection moulding, Shear viscosity, Non-equilibrium molecular dynamics, Dissipative Particle Dynamics.

**PACS:** 66.20.Cy, 05.10.Gg, 36.20.Ey

## 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 empirical laws that take into account several parameters such as 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 route 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, which arise from perturbations (linear or non-linear) of an equilibrium state, methods can be developed to treat non-equilibrium states. Furthermore, non-equilibrium molecular dynamics (NEMD) simulation methods can be devised to provide estimates for the 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

<footnote><sup>a</sup> FEMTO-ST Institute, Applied Mechanics Department, ENSMM, University of Franche-Comté, 25000 Besançon, France</footnote>

<footnote><sup>b</sup> Institute of Mechanics and Engineering, Southwest Jiaotong University (SWJTU), Sichuan, 610031 Chengdu, P.R. China</footnote>
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 reproducing hydrodynamic behaviour correctly.

## 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 20th century.

A new model, named equivalent viscosity, is proposed to take into account 4 parameters: solid volume fraction ($\phi$), shear rate ($\dot{\gamma}$), temperature ($T$) and particles size ($d_{10}$, $d_{50}$ and $d_{90}$):

$$ \frac{\eta}{\eta_0} = 10 \frac{d_{90}}{d_{10}d_{50}} A \exp(B/T) \frac{1}{\dot{\gamma}^n} \left( \frac{1-\phi}{(1-\phi/\phi_m)^2} \right) \quad (1) $$

where $\eta$ is the viscosity, $\eta_0$ the apparent viscosity, $\phi$ the solid volume fraction, $\phi_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 experimental data. The feedstock is modelled as a system composed of FeNi metallic particles and polymeric binder.

## MOLECULAR SYSTEM MODELLING

### 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 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 non-equilibrium molecular dynamics algorithm. The Green-Kubo formula for the shear viscosity is given by

$$ \eta = \frac{V}{kT} \int_{0}^{\infty} dt \langle P_{\alpha\beta}(0) \cdot P_{\alpha\beta}(t) \rangle \quad (2) $$

where $P_{\alpha\beta}$ is an off-diagonal ($\alpha \neq \beta$) of the viscous pressure tensor. $k$ is the Planck constant. $V$ and $T$ are respectively the volume and temperature of the polyatomic fluid. Angle brackets denote an equilibrium ensemble average.

$$ P_{\alpha\beta}(t) = \frac{1}{V} \left[ \sum_{i} m v_{i\alpha}(t) \cdot v_{i\beta}(t) + \sum_{i} \sum_{j>i} r_{ij\alpha}(t) \cdot F_{ij\beta}(t) \right] \quad (3) $$

where $v_{i\alpha}$ and $v_{i\beta}$ are the $\alpha$ and $\beta$ components of the velocity of particle $i$, $r_{ij\alpha}$ is the $\alpha$ component of the distance between the particles $i$ and $j$, and $F_{ij\beta}$ is the $\beta$ 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 equations of motion of the system, 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:

$$ \gamma = \lim_{X \to 0} \left[ \lim_{t \to \infty} \frac{\langle J(t) \rangle}{X} \right] \qquad (4) $$

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:

$$ \dot{\mathbf{r}}_i = \frac{\dot{\mathbf{p}}_i}{m_i} + \mathbf{r}_i \cdot \nabla \mathbf{u} \qquad (5) $$

$$ \dot{\mathbf{p}}_i = \mathbf{F}_i - \mathbf{p}_i \cdot \nabla \mathbf{u} $$

where $\mathbf{u}=(u_x,0,0)$ with $u_x = \gamma y$ is the velocity field corresponding to planar Couette flow. $\mathbf{r}_i$, $m_i$ and $\mathbf{p}_i$ are respectively the position vector, mass and momentum of particle i. $\mathbf{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, $\gamma$, cancel to yield Newton’s equations of motion relating $\mathbf{r}_i$ and $\mathbf{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 $\gamma$ 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 $\gamma$ 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

$$ \eta = \frac{P_{xy} + P_{yx}}{2\gamma} \qquad (6) $$

where $P_{xy}$ and $P_{yx}$ are the averaged xy and yx components of $\mathbf{P}$. From kinetic and mode coupling theories, to leading order the strain rate dependence of the shear viscosity is linear in $\gamma^{1/2}$. Hence, to apply the Sllod algorithm to a Newtonian fluid, one performs several simulations at differing strain rates $\gamma$ and fits the resulting strain dependent viscosities to the equation:

$$ \eta = \eta_0 + \eta_1 \gamma^{1/2} \qquad (7) $$

The zero strain rate extrapolation of $\eta$, $\eta_0$, is thus 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 facilitate 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 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 Temp in Figure 1).

<figure><img src="image_3.png" /><figcaption>FIGURE 1. Shear viscosity of feedstock composed of 60 vol% particles and 40 vol% binder at several temperatures.</figcaption></figure>

## 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 for the solving procedure of Navier-Stokes equations, and Streamline-Upwind/Petrov-Galerkin method for the solving of the transport equation [6]. Figure 2 exposes the advance of the filling front for both constant ($\mu = 10^3$ Pa.s) and equivalent viscosity formulation defined by equation (1). Global filling ratio progress is highlighted 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 ($\mu = 10^3$ Pa.s) and variable viscosity formulation (equation (1)) are highlighted. Figure 4 exposes 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.

<figure><img src="image_4.png" /><figcaption>FIGURE 2. Advance of the filling front at for significant time steps for constant (a) and equivalent viscosity (b) formulation in the case A.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIGURE 3. Evolution of the filling ratio for constant and equivalent viscosity formulations in case A.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIGURE 4. Advance of the filling front at for significant time steps for constant (a) and equivalent viscosity (b) formulation in the case B.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>FIGURE 5. Evolution of the filling ratio for constant and equivalent viscosity formulations in case B.</figcaption></figure>

CONCLUSION

NEMD is revealed a proper tool to model a particle charged fluid under shear 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.

ACKNOWLEDGMENTS

The research is supported by National Nature Science Foundation of China (No. 10772154) and New Teacher Project of Research Fund for Doctoral Program of High Education of China (No. 20070613021). The author PhD researches are funded by a scholarship supported by the Chinese Ministry of Education, the French Ministry of National Education of Higher Education and Researches (MENESR), the French Ministry of Foreign and European Affairs (MAEE) and the Conference of University Presidents (CPU). The authors would like to thank the FUI NewPIM Project.

REFERENCES

1. D. J. Evans and G. P. Morriss, "Statistical Mechanics of Non equilibrium Liquids". New York: Academic Press, 1990.

2. G. Ciccotti, D. Frenkel and I. R. McDonald, "Simulations of Liquids and Solids", 3rd ed. Amsterdam: North Holland, 1987.

3. T. Soddemann, B. Dünweg and K. Kremer, "Dissipative particle dynamics: A useful thermostat for equilibrium and nonequilibrium molecular dynamics simulations", *Physical Review E*, **68**(4), 2003.

4. P. Nikunen, M. Karttunen and I. Vattulainen, "How would you integrate the equations of motion in dissipative particle dynamics simulations?", *Computer Physics Communications*, **153**(3), pp.407-423, 2003.

5. H. J. Limbach, A. Arnold, B. A. Mann and C. Holm, "ESPResSo-an extensible simulation package for research on soft matter systems", *Computer Physics Communications*, **174**, pp.704-727 (2006).

6. G. Larsen, Z. Q. Cheng, T. Barriere, B. S. Liu and J. C. Gelin, "A Streamline-Upwind Model for Filling Front Advection in Powder Injection Moulding", NUMIFORM 2010: Proceedings of the 10th International Conference on Numerical Methods in Industrial Forming Processes Dedicated to Professor O. C. Zienkiewicz (1921-2009). AIP Conference Proceedings, **1252**, pp. 545-552 (2010).

7. G. Larsen, Z. Q. Cheng, T. Barriere, B. S. Liu, J. C. Gelin and J.C., M. R. Laydi, "Modelling and numerical simulation of biphasic fluid flow", Steel Research International, **81**(9), pp.1458-1461 (2010).

8. A. W. Lees and S. F. Edwards, "The computer study of transport processes under extreme conditions", *Journal of Physics C: Solid State Physics*, **5**(15) (1972).

9. R. Kubo, M. Toda and N. Hashitsume, "Statistical Physics II: Nonequilibrium Statistical Mechanics", 2nd ed. Springer, Berlin, 1995.