## Abstract

We present the application of the smoothed particle hydrodynamics (SPH) discretization scheme to Phillips’ model for shear-induced particle migration in concentrated suspensions. This model provides an evolution equation for the scalar mean volume fraction of idealized spherical solid particles of equal diameter which is discretized by the SPH formalism. In order to obtain a discrete evolution equation with exact conservation properties we treat in fact the occupied volume of the solid particles as the degree of freedom for the fluid particles. We present simulation results in two- and three-dimensional channel flow. The two-dimensional results serve as a verification by a comparison to analytic solutions. The three-dimensional results are used for a comparison with experimental measurements obtained from computer tomography of injection moulded ceramic microparts. We observe the best agreement of measurements with snapshots of the transient simulation for a ratio *D*_{c}/*D*_{η}=0.1 of the two model parameters.

## 1. Introduction

A great need exists for microstructured polymers, metals and ceramics. The application areas range from DNA analysis instruments in life science to biocompatible materials for medical applications to a variety of electronic devices for sensors and actuators. Injection moulding [1] represents an economically efficient process for the microstructuring of pure polymers or polymeric feedstocks with ceramic or metallic particle load. Micro powder injection moulding (MicroPIM) is a shaping process in which a highly filled polymer–powder compound (the feedstock) is heated above the melting temperature of the polymer and injected with pressure into a microcavity. The final ceramic or metallic micropart is obtained after subsequent debinding and sintering steps. To optimize part and mould design prior to tool construction and to reduce production costs there is a need for predictive process simulation.

There are a variety of finite element (FE) based simulation approaches [2,3]. Usually, their disadvantage is the difficulty in handling free surfaces and large deformations. Similarly, commercial software tools are based on one-phase material models, i.e. they show insufficient capabilities to simulate powder–binder mixtures. On the other hand, particle-based approaches are able to deal with free surfaces and large deformations easily. In this work, we apply the method of smoothed particle hydrodynamics (SPH) [4] to the simulation of the MicroPIM process. The simulation of casting with SPH was already successfully performed [5,6]. Here, we apply this method to materials which are rheologically more complex and which carry a solid load. In this work, we focus on the shear-induced segregation of the embedded powder particles [7–10] caused by the extremely high shear rates of up to 10^{6} s^{−1} in MicroPIM. This may lead to non-homogeneous binder extraction during debinding or to an anisotropic shrinkage during sintering resulting in porosity, deformations or even cracks. Barriere *et al.* [3] accounted for segregation in an FE simulation by implementing a two-phase model. This paper presents an SPH-based approach describing the powder concentration as an internal degree of freedom.

After the presentation of the governing equations for fluid flow and of Phillips’ model for shear-induced particle migration [11] in §2, we show their discretization in §3. Finally, we present several simulation results for the verification of the model, for the reproduction of experimental measurements and for complex mould geometries in §4.

## 2. Governing equations of motion

The basic equations used for the description of the MicroPIM process are the continuity equation and the incompressible Navier–Stokes equation
2.1
both given in a Lagrangian reference frame with mass density *ρ*, velocity **v**, pressure *P* and shear viscosity *η*. Segregation will be described by the diffusive flux model of Phillips *et al.* [11]. This model is based on two diffusive fluxes **J**_{c} and **J**_{η} of a mean concentration of suspended particles *ϕ*. **J**_{c} includes the migration mechanisms due to local variations in the collision frequency of the suspended particles and reads
2.2
since variations in the collision frequency are caused by concentration gradients ∇*ϕ* and variations in the shear rate . The shear rate is defined in terms of the second invariant of the strain rate tensor . *D*_{c} is a diffusion constant and *a* is the particle diameter. *D*_{c} is an empirical parameter and has to be fitted by experiments. In addition, a spatially varying viscosity owing to a spatially varying particle concentration can lead to an effective particle flux as well. It reads
2.3
where *D*_{η} is an empirical dimensionless rate constant. Using both migration mechanisms the conservation equation for the volume fraction *ϕ* becomes
2.4
The expression ∇⋅*ϕ***v** covers the convective transport of the suspended particles.

For the dependence of the viscosity on the volume fraction we assume the Krieger rheological model [12] which is
2.5
For a vanishing particle concentration the fluid has the reference viscosity *η*_{0} and for the saturation volume fraction *ϕ*=*ϕ*_{m} the viscosity diverges.

## 3. Smoothed particle hydrodynamics discretization

We briefly summarize the discretization of fluid flow. For a detailed description of the SPH discretization formalizm we refer the reader to the given references. We then show in detail how the diffusive flux model for powder migration is discretized.

To compute the evolution of the densities *ρ*_{i} of the SPH particles indexed with *i* we discretize the equation of motion in equation (2.1), i.e.
3.1
where we have set **v**_{ij}=**v**_{i}−**v**_{j} and ∇*W*_{ij}=−**r**_{ij}*w*_{ij}. *W*_{ij} is the SPH interpolation kernel and *m*_{i} a particle mass. This equation allows one to simulate free surface flow [13]. The discretization of the pressure gradient term in equation (2.1) reads
3.2
Incompressibility is approximated by a weakly compressible equation of state [13]
3.3
where *γ*=7, *ρ*_{0} is a reference density and *P*_{0}=*c*^{2}*ρ*_{0}/*γ*, with the speed of sound *c* chosen large enough so that density fluctuations remain below 1 per cent. The viscous term in the momentum equation in (2.1) is discretized as [14]
3.4
Finally, the tensor ∇** v** can be discretized by [4]
3.5
which we require for the computation of the shear rate in equations (2.2) and (2.3).

Instead of discretizing equation (2.4) directly, it is better to use the equation of motion for the volume *V*_{ϕ}=*ϕ*/*ρ* occupied by powder particles, where we have set the mass to *m*=1*m**, i.e. to the unit of mass *m** for simplicity. The discretized equation for *V*_{ϕ} will be antisymmetric under exchange of fluid particles, leading to exact volume conservation. We have to switch from the Eulerian to the Lagrangian description. Introducing the material derivative d/d*t*=∂/∂*t*+**v**⋅∇ and substituting *ϕ*=*ρV*_{ϕ} gives
3.6
which is the equation of motion for the occupied volume. The divergence of the velocity field gives the relative volume expansion, i.e. . Thus, the first two terms on the right-hand side of (3.6) cancel. Now, we apply the rule for the SPH discretization of second derivatives [4]
3.7
for two scalar fields *A*(**r**) and *B*(**r**) to the right-hand side of equation (3.6). This is a simpler form than the discretization (3.4), but we have not found any significant difference between the two forms for the results shown here. We get
3.8
Note that the expression within the sum is antisymmetric under particle exchange. In the discretized equation for , a factor 1/*ρ*_{i} is missing and the equation has no symmetry. The advantage of the antisymmetry is that volume is exactly conserved, because the amount of suspended volume that is deducted from particle *i* is added to particle *j*.

We finally have to consider the numerical implementation of the Krieger model. Both *η*(*ϕ*) and the derivative (d*η*/d*ϕ*)/*η* from equations (2.5) diverge for *ϕ*→*ϕ*_{m}. Therefore we have to use approximate expressions giving finite values for concentrations *ϕ*_{c}<*ϕ*<*ϕ*_{m}, where *ϕ*_{c} is a critical limit. It is most convenient to start with an approximation for the relative derivative (d*η*/d*ϕ*)/*η*. The lowest order approximation we can make is
3.9
where *η*_{d}(*ϕ*) represents the form of the Krieger model used for the discretized equations and may be obtained by solving the differential equation given in equation (3.9), i.e.
3.10
For a continuous transition to *η*(*ϕ*) at *ϕ*=*ϕ*_{c} we need the boundary condition *η*_{d}(*ϕ*_{c})=*η*(*ϕ*_{c}). Then, solving equation (3.10) gives
3.11
For *ϕ*=*ϕ*_{m} a finite value *η*_{d}(*ϕ*_{m})=*η*(*ϕ*_{c})e^{c} is obtained. Furthermore, since the approximation for the derivative (3.9) does not diverge for *ϕ*→*ϕ*_{m} the flux **J**_{η} from equation (2.3) is limited artificially and allows for particle concentrations *ϕ*>*ϕ*_{m} above the packing limit. An unphysical effect is to be expected only when and *ϕ*>*ϕ*_{c}. For the steady-state simulations presented in this work it turns out that the smoothing effect of the SPH discretization (3.5) of counteracts and dominates over this approximation. In the transient simulations—which resemble the PIM process more closely—the case *ϕ*>*ϕ*_{c} never occurs.

## 4. Simulation results and discussion

The SPH form of the Phillips model is verified against an analytic solution by computing the concentration profile for pressure-driven Poiseuille flow between two parallel plates. We choose *d*_{x}=250 μm as the width of the cross section in *x*-direction. This width corresponds to the experimental set-up for which measurements from computer tomography (CT) are available. The velocity component *v*_{y} is driven by a constant pressure gradient *p*_{y}=1.28 MPa cm^{−1} which is implemented by applying an equivalent body-force on each SPH particle in the periodic flow direction. The particle diameter is taken as *a*=0.01 μm, *D*_{c}=1.17×10^{9} is chosen such that the time constant for reaching steady state is as small as possible and we set *D*_{c}/*D*_{η}=(1+1/*c*)^{−1}≈0.645. This allows for an analytic solution of the steady-state concentration profile which is for the left half of the channel *x*^{′}=[0,1/2]
4.1
with *x*^{′}≡*x*/*d*_{x} and where *ϕ*_{0} is the concentration at the walls. Here, we take *ϕ*_{0}=0.45, *ϕ*_{m}=0.68 and *c*=1.82 in accordance with Phillips *et al*. [11].

The simulation is performed in three dimensions with periodic boundary conditions in the flow direction *y* and in *z*-direction. The cross section in *x*-direction is discretized with two different resolutions of 9 and 19 SPH particles which are initially at rest on a cubic lattice. The *x*-direction is confined by solid boundaries where no-slip boundary conditions are applied by boundary particles using the technique of Morris *et al*. [15]. By forbidding any exchange of *ϕ* between the SPH particles representing the liquid and the SPH particles representing the solid walls we ensure zero migration across the boundaries. Figure 1 shows the powder concentration, the velocity and the shear rate profiles in steady state and the analytic solution (4.1). One observes that the solid migrates to the centre of the channel and a blunted velocity profile. The latter is due to an increased viscosity in the centre because of the larger solids concentration. Hence, it is crucial to add the degree of freedom *ϕ* to account for its effect on the rheology. The computed concentration does not fully reach the analytic result *ϕ*=*ϕ*_{m}=0.68 at the peak owing to the smoothing of SPH. Thus, mass conservation leads to concentrations at the walls larger than *ϕ*_{0}=0.45 when compared with the analytic solution. The error decreases when increasing the resolution.

Figure 2 shows the transient evolution of the powder concentration as obtained from a simulation of flow through a quadratic 250×250 μm channel cross section. The simulation settings are the same as before except for *D*_{c}/*D*_{η}=0.4, *a*=9 μm, *η*_{0}=33.2 Pa s, *ϕ*_{m}=0.794 and *c*=0.714, which are taken now from viscosity measurements of our feedstock ‘GoMicro’ [9]. This corresponds to a ‘mean’ viscosity , where . We have simulated different ratios *D*_{c}/*D*_{η} and have found that *D*_{c}/*D*_{η}=0.1 reproduces best our CT measurement for a specially prepared injection-moulded test specimen [7,9]. For this case, figure 3*a* plots the concentrations at the centre line of the cross section and also shows a concentration profile as obtained from CT and, for comparison, a profile obtained from a simulation with *D*_{c}/*D*_{η}=0.2, showing the trend when increasing *D*_{c}/*D*_{η}. Figure 3*b* shows the CT data in two dimensions. The global average concentration was conserved exactly in the simulations. In the simulations presented here, the equilibration of *ϕ* was boosted on purpose by choosing large constants *D*_{c} and *D*_{η}. Segregation starts close to the walls and propagates towards the centre where, after *t*=0.1 *t** a cusp starts to grow. Around this time the profile is similar to our CT measurement. This indicates that steady state was not reached during the MicroPIM process, i.e. the longer the channel, the larger the segregation. Observe also the similarity between figures 2*b* and 3*b*. It is known that the ratio of diffusion constants *D*_{c}/*D*_{η} is material dependent. Our value of *D*_{c}/*D*_{η}=0.1 is small with respect to the physically reasonable range [0..1] [11]. Hence, the flux **J**_{η} from equation (2.3) is dominant, which is a reasonable observation for highly viscous PIM feedstocks.

As a more complex application of the powder migration model we consider the simulation of injection moulding into a three-dimensional geometry. Some of our real PIM parts are shown in figure 4*a* and a snapshot of the filling simulation in figure 4*b*. The arrows in figure 4*b,d* indicate the flow direction at the gate. Besides the similarity of the flow patterns, figure 4*a,b* illustrates the advantages of SPH over FE methods. The SPH particles naturally reproduce the large deformations and the free-surface-dynamics of the feedstock material which do not have to be extracted from a fixed grid. By the SPH discretization of Phillips’ model presented in this work, it is now possible to track directly the solid load for which no CT data are available. This is shown in figure 4*c*,*d* which focuses on the gate region. Figure 4*c* is a snapshot of an early stage of the injection. The distribution of the solid load indicates an aggregation at convex corners (i.e. pointing into the mould material) and a decrease at concave corners. This is intuitively understandable since concave corners are regions of large shear rates. In later stages, the flow close to the gate becomes more homogeneous and directed towards the two arms at the top and the bottom. In these narrow arms strong shear rates occur. The effect can be seen in figure 4*d*. The concentration in the entrance volume is rather homogeneous and large. Inside the arms the concentration is lower. A larger fraction of the solid particles prefers to stay in the entrance volume where shear rates and flow velocities are lower. Additionally, the solid load has a maximum in the centre of the arms as already observed earlier for flow in a channel.

## 5. Summary

An SPH framework was developed which allows one to simulate the injection process of MicroPIM. The simulations can be performed in arbitrarily complex three-dimensional geometries. Shear-induced powder migration was incorporated by means of Phillips’ diffusive flux model. This model was discretized by formulating an SPH equation of motion for the occupied volume *V*_{ϕ} with exact pairwise conservation properties. The simulations correctly predict powder migration to regions with the lowest shear rates and indicate that, for the available CT measurements, the measured profile is not a steady state. For injection moulding into complex geometries the simulations predict regions of accumulation of the solids fraction. For a quantitative matching of the transport coefficients *D*_{c} and *D*_{η} of the Phillips model we suggest CT measurements for well defined geometries such as long capillaries.

## Acknowledgements

The authors acknowledge funding of this project by the Deutsche Forschungsgemeinschaft (DFG) in the framework of the Sonderforschungsbereich Mikrourformen (SFB499) and by the University of Freiburg through the German excellence initiative.

## Footnotes

One contribution of 26 to a Theme Issue ‘Discrete simulation of fluid dynamics: methods’.

- This journal is © 2011 The Royal Society