Analysis of Powder Compaction Process Using Multi-Particle Finite Element Method

Abstract

Cold compaction is an essential step of the powder metallurgy process, which is cost-effective in the manufacturing industry. In the present study, the compaction of pure iron powder is investigated. A MATLAB code is developed to create a representative volume element with a given size distribution, number of particles, and initial relative density. The finite element analysis is performed by implementing an ABAQUS/Explicit python script. Multi-particle finite element analysis of compaction is employed to analyze the deformation of the particles into the green body. The effects of loading path, geometry configuration, and wall friction on punch force, axial and transverse stresses, strains, and yield surface are discussed. The results indicate that wall friction affects the load, stress state, and yield surface. It was found that, with the increase of the wall friction from 0.0 to 0.2, the compaction force increases by 14.70%. Also, the difference between the upper and lower RVE face forces increases from 1.01% to 16.27%. Changing the loading path from the compression to the hydrostatic compaction decreases the compaction force by 57.7% in upper and 65.75% in lower RVE faces, while the axial stress decreases by 27.70% and the transverse stress increases by 40.41%. The volume strain of compaction is smaller by 35.33% than that of hydrostatic compression.

1. Introduction

Powder metallurgy is helpful for the mass production of parts with complex geometries. An essential step of this process, cold compaction of powders, is of interest for many manufacturing processes, for example, the production of green bodies, ceramic parts, pellets, and tablets in the pharmaceutical industry.^1,2) In lots of uses, powder parts must be compacted up to as high as possible relative density.³⁾

The finite element method is a powerful tool to model the process and predict the process outputs. Many researches used numerical methods to investigate the flow behavior in the powder compaction process. Three main modeling approaches in powder compaction simulations are the continuum media approach, the discrete element method (DEM), and the multi-particle finite element method (MPFEM).

Inter-particle contact interactions were ignored in the continuum media approach, and particles were studied as a continuous medium. In this method, the macroscopic behavior of the material can be investigated as a plasticity problem with significant deformation.⁴⁾ Employing this approach, Sinka et al. investigated the effect of wall friction in pharmaceutical tablets compaction using the Drucker–Prager cap model.⁵⁾ Redanz et al. and Lewis et al. used the continuum approach to analyze large deformation in the compaction process.^6,7) Akhmetshin et al. simulated compaction of tungsten carbide/cobalt powder via the Drucker-Prager Cap constitutive model to obtain the stress distribution and volume of plastic deformation during compaction.⁸⁾

Cundall⁹⁾ developed DEM in 1971. In this method, the motion of particles is studied through the Newton’s second law of motion, considering contact interactions between particles and processing contact forces as functions of the particle overlap.^10,11) In this regard, Jeong and Choi used DEM to investigate the stress transmission and compressive force distribution between the top and bottom punches in metal powder compaction with the addition of graphite powders.¹²⁾ Zhang et al. utilized DEM to investigate the correlation mechanism between force chains and the friction mechanism during powder compaction.¹³⁾ Xu and Meng employed DEM to study metal powder flow and force chain characteristics in high-velocity compaction.¹⁴⁾ Redanz and Fleck used DEM to investigate the yield surfaces for sticking and frictionless particle contacts in a 2-dimensional analysis.¹⁵⁾

In MPFEM, classical finite element simulations on a discrete particle assembly are performed. These particles are composed of finite elements, and their interactions are ruled by classical finite-element contact conditions.¹⁶⁾ Using MPFEM, Abdelmoula et al. investigated the plastic flow of granular materials with highly deformable elastic-plastic grains.¹⁷⁾ Demirtas and Klinzing investigated the compaction behavior of hollow spheres computationally and experimentally.¹⁸⁾ They simulated two different particle sizes with different shell thicknesses using MPFEM.

The most realistic MPFEM simulation that considers all the affecting parameters has been done with a limited number of particles so far. This research aims to investigate the macrostructure of the powder compaction process using a vast number of particles. In section 2, the representative volume element generation method is explained, followed by the finite element method. Section 3 describes the results and discussion. Section 4 concludes the research achievements.

2. Methods

2.1 Material properties

In this research, the elasto-plastic model has been used to analyze pure iron powder compaction. Table 1 includes the material properties for the elasto-plastic model. The stress-strain equation for this model is as follow:   

\begin{equation} \bar{\sigma} = 465.5(\bar{\varepsilon} + 0.0140)^{0.2481} \end{equation} (1)

where $\bar{\sigma }$ (in MPa) and $\bar{\varepsilon }$ are the effective stress and the equivalent plastic strain, respectively. The true stress-strain diagram for both material models is shown in Fig. 1.

Table 1 Material parameters of pure iron powder.

Fig. 1

Stress-strain curve of the iron powder.

2.2 Representative volume element generation

The first step of the simulation is to generate a representative volume element (RVE) that consists of a desired number of particles considering initial relative density and preventing overclosure between particles. A MATLAB code has been implemented to generate the RVE. Table 2 includes the inputs and outputs of the RVE generator code.

Table 2 Input parameters of the RVE generator.

Figure 2 shows the experimentally measured particle size distribution and the output distribution. The curve had been fitted using a Gaussian distribution considering the cut-off radii (22.5 µm∼88.5 µm):   

\begin{equation} y = \frac{1}{2\pi\sigma}e^{-\frac{(x - \mu)^{2}}{2\sigma^{2}}} \end{equation} (2)

where μ is the mean of the distribution and σ is the standard deviation, taking the values of 111 µm and 29 µm in this work to fit the experimental distribution, respectively.

Fig. 2

Measured and generated particle size distribution.

The algorithm for generating RVE is as follows:

1. 1.    Determine the RVE size based on the Gaussian fitted particle size distribution and input parameters.
2. 2.    Determine the RVE edge size, a:   
   \begin{equation} a = \left[\frac{4\pi r_{m}^{3}}{3\rho_{i,\textit{target}}}\right]^{1/3} \end{equation} (3)

where ρ_i,target is the target initial relative density, and r_m is the mean radius.

1. 3.    Generate a list of radii with cut-off radii.
2. 4.    Choose a radius.
3. 5.    Divide the RVE into G³ cubes of an equal size (G: an integer number).
4. 6.    Sweep all the cubes in a random order and check if a new particle space is occupied by other particles. For each cube, maximum 50 new positions are tried. If an unoccupied space is found, insert the particle and go to Step 8.
5. 7.    If the sweeping is completed but fails to find a new particle position, increase G by 2 and go to Step 5.
6. 8.    Repeat 4–7 until no more particles can be inserted.

2.3 Finite element model

Using the generated RVEs, finite element models were constructed for investigating the effect of various model and process parameters. Table 3 shows the outputs of the generated RVE for FEM studies.

Table 3 Output parameters of the RVE generator.

The finite element modeling and analysis have been done by ABAQUS/Explicit 2020 python script. Figure 3 shows a generated RVE.

Fig. 3

(a) The RVE with loading directions, and (b) the discretized RVE.

2.3.1 Model sensitivity analysis

To check the sensitivity of the finite element mesh, the effect of model parameters has been examined. Figure 4(a) and (b) show the effect of the number of elements and element type on the maximum punch force and analysis time for the elasto-plastic model, using linear and quadratic tetrahedral elements.

Fig. 4

Sensitivity analysis of the FE model parameters: (a) the effect of element number on the analysis time, (b) the effect of element number on the maximum punch force, (c) the effect of mass scaling factor on the pressure-deviatoric stress (p-q) diagram, and (d) the effect of punch speed on the average equivalent plastic strain rate.

Figure 4(c) shows the p-q diagram for different mass scaling factor. The pressure (hydrostatic stress) (p) and the equivalent stress (q) are computed as functions of the main principal stresses of the stress tensor at different compaction stages:   

\begin{equation} p = \frac{1}{3}(\sigma_{11} + \sigma_{22} + \sigma_{33}) \end{equation} (4)

  

\begin{equation} q = \sqrt{\frac{1}{2}[(\sigma_{11} - \sigma_{22})^{2} + (\sigma_{33} - \sigma_{22})^{2} + (\sigma_{11} - \sigma_{33})^{2}]} \end{equation} (5)

where σ₁₁, σ₂₂, and σ₃₃ are the three normal stress components.

Figure 4(d) shows the average equivalent plastic strain rate ($\dot{\bar{\varepsilon }}$ or ER_eq) versus punch speed curve. ER_eq can be calculated as follow:   

\begin{equation} \dot{\bar{\varepsilon}} = \frac{2}{3}\sqrt{\frac{1}{2}[(\dot{\varepsilon}_{1} - \dot{\varepsilon}_{2})^{2} + (\dot{\varepsilon}_{2} - \dot{\varepsilon}_{3})^{2} + (\dot{\varepsilon}_{3} - \dot{\varepsilon}_{1})^{2}]} \end{equation} (6)

where $\dot{\varepsilon }_{1},$ $\dot{\varepsilon }_{2} $, and $\dot{\varepsilon }_{3} $ are the three principal strain rates.

Considering all the described effects, Table 4 displays the selected configurations for the finite element model.

Table 4 Selected configurations of the finite element model.

3. Results and Discussion

3.1 Effect of initial configuration

Ten different initial configurations (based on experimental size distribution) have been analyzed to study the influence of geometry. Table 5 shows the properties of generated initial configurations. Figure 5 shows the ten different RVEs. The edge length of the RVEs is 1015.9 µm.

Table 5 Information of the initial RVE configurations.

Fig. 5

Initial RVE configurations: (a), (b), …, (j) correspond to RVE1, RVE2, …, RVE10, respectively.

Figure 6 shows the effect of different initial configurations on the punch force, p-q diagram, and the average logarithmic strain distribution frequency in the compaction, while Fig. 7 shows the effect of the initial configuration on the average logarithmic strain distribution frequency.

Fig. 6

Effect of different initial RVE configurations on (a) the punch force, and (b) the p-q diagram.

Fig. 7

Average logarithmic strain for different initial RVE configurations.

As can be seen, the effect of different initial configurations on punch force and is negligible. The standard deviation of the maximum punch force equals 5.22 kN, which is about 2.2% of the average maximum punch force (229.73 kN). On the other hand, some deviations can be seen in the p-q diagram, which shows the initial configuration affects the stress state.

3.2 Effect of wall friction

Three different wall friction has been analyzed in the model to investigate the effect. Figure 8(a) shows the upper and lower RVE face force of compaction die in different wall frictions. As can be seen, the wall friction affects the punch force, With the increase of the wall friction from 0.0 to 0.2, the compaction force increases by 14.70%. Also, the difference between the upper and lower RVE face forces increases from 1.01% to 16.27% when the wall friction coefficient rises.

Fig. 8

Effect of wall friction on: (a) the punch force, (b) the stress components, and (c) the p-q diagram.

Figure 8(b) shows the stress components of compaction RVE in different wall frictions. Axial and transverse stresses also slightly increase with the increase of friction coefficient, which is due to the rise of the forces because of friction on the stresses.

Figure 8(c) shows the p-q diagram for these cases. When the relative density is about 0.33%, the q in the p-q diagram is slightly smaller when the friction coefficient equals zero, but this difference is negligible when the relative density increases to 72%.

3.3 Effect of loading path

Twenty-one different speed ratios have been simulated to investigate the loading path effect considering the loading path parameter α (the ratio of the transverse speed to the axial speed, Fig. 3) equal as −1.0, −0.9, −0.8, …, 0.0, 0.1, 0.2, …, 1.0 (Fig. 3). Figure 9(a) shows the punch force versus relative density diagram for three loading paths (α = 0, 0.5, 1.0), respectively.

Fig. 9

Effect of loading paths on the punch force: (a) the punch force, and (b) the stress components.

When α increases from 0 to 1, the loading path changes from compaction to hydrostatic compression. It can be seen that punch force decreases by 57.70% in upper and 65.75% in lower RVE faces with the increase of the α. The loading path with α = 0 has the highest punch force and axial stress.

Figure 9(b) shows the axial and transverse stresses versus the relative density diagram of three loading paths (α = 0, 0.5, 1.0), respectively. When the α changes from 0 to 1, the decrease of axial stress (27.70%) accompanied by an increase in the transverse stress (40.41%) leads to the increase of the hydrostatic stress.

To better understand the morphology of particles, the initial and final ZX section views of particles are shown in Fig. 10, along with the void percentage and average volume strain. The percentage of void marginally decreases with the increase of the α. It can be concluded that the loading path affects the plastic deformation; with the increase of α, the amount of plastic deformation increases.

Fig. 10

Initial and final ZX section views of particles and average volume strains at the relative density of 71% for three loading paths: α = 0, 0.5, and 1.0.

To investigate the springback phenomena in the compaction process, the average logarithmic strain of particles is calculated in three different load paths: α = 0.0, α = 0.5, and α = 1.0 (Fig. 10). The average volume strain by 35.33% with the increase of α.

As shown in Figs. 11(a) and (b), the p-q diagram and the yield surface change dramatically with the change in loading path; the main difference exists when α = 0. The average of three different initial configurations has been considered to calculate the yield surface. Figure 11(c) is a histogram of the average logarithmic strain that shows less strain frequency with increasing α. Also, it can be seen in Fig. 11(d) that level of logarithmic strain decreases with increasing α.

Fig. 11

(a) The p-q diagram, (b) the yield surface, (c) the average logarithmic strain distribution, and (d) the average logarithmic strain of individual particles for different loading paths.

4. Conclusion

Multi-particle finite element simulations of the powder compaction process have been performed. Modules of particle generation, FE mesh generation, FE calculation, postprocessing were developed using MATLAB, Python, and ABAQUS. An initial relative density of 26.7% was achieved with 456 particles. The effect of model and process parameters has been examined. The following points can be concluded:

1. (1)    Linear tetrahedral elements with similar sizes to the mean radius were found to be effective in MPFEM.
2. (2)    The effect of different initial configurations on punch force is negligible. Initial configuration affects the stress state. When relative density is about 40% of the p-q diagram, the q is slightly smaller when the friction coefficient is zero. Still, this difference is negligible when the relative density increases to the maximum value.
3. (3)    With the increase of the wall friction from 0.0 to 0.2, the compaction force increases by 14.70%.
4. (4)    The difference between the upper and lower RVE face forces increases from 1.01% to 16.27% when the wall friction coefficient increases from 0.0 to 0.2.
5. (5)    When the loading path changes from compaction to hydrostatic compression, the punch force decreases by 57.70% in upper and 65.75% in lower RVE faces.
6. (6)    When the loading path parameter or the velocity ratio, α, increased, the decrease of axial stress (27.70%) accompanied by an increase in transverse stress (40.41%), leads to the increase of the hydrostatic stress.
7. (7)    The average logarithmic strain decreases with the increase of α.
8. (8)    The percentage of void marginally decreases from 28.48% to 28% with the increase of the α due to increased plastic deformation.
9. (9)    The volume strain of compaction is more significant than that of hydrostatic compression and decreases by 35.33% with the increase of α.
10. (10)    The yield surface was successfully calculated utilizing different loading paths.

Acknowledgments

This work was supported by the National Research Foundation of Korea (NRF) Grant funded by the Korean Government (MSIP) (No. Grant Number – NRF-2021M3H4A6A01045764).

REFERENCES

• 1)   S.C.  Lee and  K.T.  Kim: Int. J. Mech. Sci. 44 (2002) 1295–1308. doi:10.1016/S0020-7403(02)00054-1
• 2)   M.  Trunec and  K.  Maca: J. Am. Ceram. Soc. 90 (2007) 2735–2740. doi:10.1111/j.1551-2916.2007.01781.x
• 3)   P.  Loidolt,  M.H.  Ulz and  J.  Khinast: Powder Technol. 336 (2018) 426–440. doi:10.1016/j.powtec.2018.06.018
• 4)   I.C.  Sinka: Kona Powder Particle J. 25 (2007) 4–22. doi:10.14356/kona.2007005
• 5)   I.C.  Sinka,  J.C.  Cunningham and  A.  Zavaliangos: Powder Technol. 133 (2003) 33–43. doi:10.1016/S0032-5910(03)00094-9
• 6)   P.  Redanz: J. Mech. Solids 18 (1999) 399–413. doi:10.1016/S0997-7538(99)00120-5
• 7)   R.W.  Lewis and  A.  Khoei: Comput. Methods Appl. Mech. Eng. 159 (1998) 291–328. doi:10.1016/S0045-7825(97)00276-4
• 8)   L.R.  Akhmetshin,  I.Y.  Smolin and  S.P.  Buyakova: Mater. Sci. Eng. 1093 (2021) 012002. doi:10.1088/1757-899X/1093/1/012002
• 9)  P.A. Cundall: Proceedings of the Symposium of the International Society for Rock Mechanics, Nancy, France, 1, (1971) pp. 2–8.
• 10)   G.  Gustafsson,  H.-Å.  Häggblad and  P.  Jonsén: Powder Technol. 239 (2013) 231–238. doi:10.1016/j.powtec.2013.02.005
• 11)   P.  Cundall and  O.  Strack: Geotechnique 29 (1979) 47–65. doi:10.1680/geot.1979.29.1.47
• 12)   J.H.  Jeong and  J.  Choi: Powder Metall. 64 (2021) 241–247. doi:10.1080/00325899.2021.1921404
• 13)   N.  Zhang,  Sh.  Zhang,  J.  Tan and  W.  Zhang: Chin. Phys. B 31 (2022) 024501. doi:10.1088/1674-1056/ac0db1
• 14)   P.  Redanz and  N.A.  Fleck: Acta Mater. 49 (2001) 4325–4335. doi:10.1016/S1359-6454(01)00298-1
• 15)   Z.Y.  Xu and  F.J.  Men: J. Korean Phys. Soc. 79 (2021) 455–467. doi:10.1007/s40042-021-00241-9
• 16)   B.  Harthong,  J.F.  Jérier,  V.  Richefeu,  B.  Chareyre,  P.  Dorémus,  D.  Imbault and  F.V.  Donzé: Int. J. Mech. Sci. 61 (2012) 32–43. doi:10.1016/j.ijmecsci.2012.04.013
• 17)  N. Abdelmoula, B. Harthong, D. Imbault and P. Doremus: IV International Conference on Particle-based Methods–Fundamentals and Applications, (2015).
• 18)   A.  Demirtas and  G.R.  Klinzing: Powder Technol. 391 (2021) 34–45. doi:10.1016/j.powtec.2021.06.004