# A Study of Compaction Densification Behavior of Composite Particles by Multiparticle Finite Element Method

Authors: Junnan Li; Peng Han; Yingshi Xu; Kun Liu; Guoxin Zhang; Peng Xu
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

In this paper, 3D particulate scale modelling on the die compaction of DEM generated initial packings of both soft and hard particles was conducted by employing the multiparticle finite element method (MPFEM). The effects of initial packing structures as well as the compaction pressure on the macroscopic and microscopic properties of the whole powder mass and local structures were investigated. In addition, corresponding physical experiments were carried out for model validation. The results show that the compact obtained from the initial dense packing under vibration undergoes yielding stage earlier than that with natural initial packing (without vibration), and the relative density is much higher. Pores that are significantly smaller and with more uniform size and homogenous stress distribution are observed in the former case. Highest stress regions occur in most cases at a grain boundary with large curvature after deformation. Moreover, the high stress in the central part of both soft and hard particles during compaction is significantly reduced after pressure unloading, reaching a new force balance. In this case, the stress is concentrated mainly at the corners of the deformed particles, which creates the risk of cracking during subsequent sintering at either the contact region between particles or the corners. The numerical results are found to be in good agreement with those from physical experiments, confirming the robustness and reliability of the numerical model used in the simulations.


---

## Body

## 1. Introduction

Powder metallurgy (PM) has long been an indispensable technology for both advanced and conventional materials [1–4]. Its advantages lie in cost-effectiveness, straightforward operation, conservation of materials, environmental friendliness, high production rates, and net or near net shape forming for complex geometrical and structural metals or alloys with high melting point or porous structure. These properties make it very competitive with conventional metallurgical processes such as casting or other processes such as forging, machining, and others [5, 6]. An important part of the preparation of PM for these materials is powder forming, which has a large impact on the subsequent sintering process and on the performance of the final PM product. For this reason, a great deal of research attention, both experimental and theoretical, has been focused on this in the previous decades.Regarding physical experiments, most of the focus is on establishing the relationship between the relative density of the compaction and the external pressure during the powder forming process. Since Walker first proposed the linear relationship between the relative volume of the powder and the logarithm of the compaction pressure based on experiments, various formulas have been proposed to predict the quantitative relationship between compaction pressure and relative density. However, it proved difficult to quantitatively characterize the local density distribution, stress distribution, and particle flow evolution process of the green body in the powder forming process with physical experiments, due to the geometric, material, and contact nonlinearities [7–11]. In this context, various numerical models describing the forming of metal powders have been proposed. From the perspective of metal plastic mechanics, based on the von Mises theory of metal materials, Kuhn [12], Green [13], Shima [14], Gurson [15], Doraivelu [16], and others have proposed different forms of yield criteria for powder materials. For soil plastic mechanics methods, Drucker-Prager-Cap model [17] and Cam-Clay model [18] are mainly represented. Based on these yield models, a macrocontinuous finite element method (FEM) was developed to solve the difficult problems encountered in physical experiments. The main advantage of the macroscopic continuous scale finite element method is that it can be used to analyze the macroscopic characteristics of the local density and distribution, the stress distribution, and the particle displacement during the forming process. In previous research, this method was successfully used for unidirectional molding of pure metal powder [19–21] and composite powder [22, 23]. Although the traditional finite element method can avoid the difficulties present in physical experiments to a certain extent, the powder particles are regarded as a continuum with uniform distribution in the FEM model, so it is difficult to characterize the particle dynamics and contact mechanical behavior in the powder forming process from the microscale. The deficiency of FEM simulation in this aspect can be compensated with a discrete element (DEM) simulation. Using this method, An et al. [24–27] successfully simulate the packing densification of spherical and nonspherical particles. However, the study found that the limitation of DEM is that it is only suitable for the simulation of small deformation particle accumulation, and it is difficult to simulate the formation of particles with large plastic deformation; in particular when the packing density of particles exceeds 0.85, it is very difficult to implement DEM [28]. Moreover, since it can amend the drawbacks presented by both the FEM and DEM simulation, the multiparticle finite element method (MPFEM) proposed in recent years has been gaining currency [29–43]. However, few studies on the simulation of soft and hard particles forming process from the particle scale using multiparticle finite element (MPFEM) can be found in literature. Only Wu et al. [40] simulated the 2D forming process of composite powder, and only the formation of ordered initial packing structure was considered. In fact, the initial packing of powder before forming is usually disordered, and the research in this aspect is only considered in our previous 2D simulation of composite powder [44, 45]. However, the simulation of 2D forming process is not sufficient to describe the evolution of compact structure and properties as well as the micromechanism of densification during 3D powder forming, and the research is not presently available.For this reason, in this paper, MPFEM was used to simulate the compaction process of soft and hard particles on a 3D particle scale. The influence of powder density and initial packing structure was considered throughout the entire process. The density and distribution of compacts, stress and distribution, morphology evolution of particles, filling mechanism of pores, stress transfer behavior and variation of internal stress of different particles during forming, and pressure relief of soft and hard particles were systematically analyzed. The results are found to be of great significance for the preparation of high-performance soft and hard particles compacts and the analysis of the densification mechanism.


## 2. Simulation Method and Conditions

### 2.1. Model Setup

The particle scale MPFEM numerical simulation was carried out using commercial MSC Marc software. In order to reduce constraints and obtain a better simulation of the physical process, DEM dynamic modelling was employed to generate initial random packing structure which was then imported into the FEM model for compaction. The DEM method has been successfully applied to packing of various spherical and nonspherical particulate systems [24–27]. Figure 1(a) shows the DEM generation of initial random packing structure with discretized mesh division on all particles before compaction when the mass fraction of soft particles is Xsoft = 25wt.%, where each particle includes 1014 units and 872 nodes as shown in Figure 1(b). In this case, both soft and hard particles have equal size (1 mm in diameter). At the interface between two particles in contact, the Coulomb friction model was used:(1)σfr=−μσn2πarctanvrrvcnstt,where sfr is the tangential (friction) stress; sn is the normal stress; µ is the friction coefficient; t is the tangential vector in the direction of the relative velocity; vr represents the relative sliding. Physically, the value of rvcnst can be regarded as the value of the relative velocity below which sticking occurs.Figure 1
DEM generated 3D initial random dense packing structure (a) and mesh division of an individual particle in the soft and hard particles (b), the packing structures for with vibration (c) and (d) without vibration.
(a)(b)(c)(d)In the simulation, the yield stress can be expressed by(2)σy=Cε0+ε¯n+Dε¯˙k,where e0 is initial yield strain; ε¯ is equivalent strain; ε¯˙ is equivalent strain rate; parameters C, D, n, and k are material constants. Initially, the equivalent strain and equivalent strain rate were assumed to be zero so that the initial yield strain can be calculated by material constants C and n and Young’s modulus Y in von Mises material model. The von Mises stress is given by(3)σ¯=σ1−σ22+σ2−σ32+σ3−σ121/22,where s1, s2, and s3 are the principal Cauchy stresses along three main axes. The equivalent strain rate is(4)ε¯˙=23εijεij1/2,where eij is the strain tensor.MSC Marc uses Cauchy stress (true stress) and logarithmic strain with updated Lagrange formulation. It is instructive to derive the stiffness matrices for the updated Lagrangian formulation starting from the virtual work principle in the following equation:(5)∫V0SijδEijdV=∫V0bi0δηidV+∫A0ti0δηidA.Here,Sij is the second Piola-Kirchhoff stress; Eij is the Green-Lagrange strain; bi0 is the body force in the reference configuration; ti0 is the traction vector in the reference configuration; and ηi is the virtual displacement. Direct linearization of the left-hand side of (5) yields(6)∫V0SijdδEijdV=∫Vn+1∇ηikσkj∇Δuijdv,where Δu and η are actual incremental and virtual displacements, respectively, and skj is the Cauchy stress tensor. Equation (6) can be further calculated as(7)∫V0dSijδEijdV=∫Vn+1∇sηijLijk1∇sΔuk1dv,where ∇S denotes the symmetric part of ∇, which represents the gradient operator in the present configuration. Lijk1 is convected to the present configuration to yield:(8)K1+K2δu=F−R,where F and R are the external and internal forces, respectively; δu is the correction displacement vector; K1 and K2 are, respectively, the material stiffness matrix and geometric stiffness matrix, given by(9)K1ij=∫Vn+1βimnLmnpqβpqjdv,and(10)K2ij=∫Vn+1σk1Ni,kNj,1dv,where ßimn is the symmetric gradient operator evaluated in the present configuration; Lmnpq represents the material moduli tensor; sk1 is the Cauchy stresses; Ni,k is the shape function gradient matrix.Keeping in view that the reference state is the present state, a rate formulation can be expressed as(11)∫V0S˙ijδEij+Sij∂vk∂Xi∂δηk∂XjdV=∫V0b˙iδηidV+∫A0t˙iδηidA,where S˙ij is the rate of the second Piola-Kirchhoff stress. Equation (11) can be obtained by setting(12)Fij=δij,δEij=δdij,∂∂Xi=∂∂xi,Sij=σij,where Fij is the deformation tensor; dij is the rate of deformation; Xi is the coordinate of initial node; and xi is the coordinate of the node. Hence:(13)∫Vn+1σ∇ijδdij+σij∂vk∂xi∂δηk∂xjdv=∫Vn+1b˙iδηidv+∫Vn+1t˙iδηida,in which b˙i and t˙i are the body force and surface traction in the present configuration, respectively. And in this equation, σij∇ is the Truesdell rate of Cauchy stress.MSC Marc uses Cauchy stress (true stress) and logarithmic strain with updated Lagrange formulation (explained in more detail in the preliminary work [44, 45]).


### 2.2. Simulation Conditions

The boundary conditions of the MPFEM simulation are given as follows: displacement boundary condition is imposed via rigid upper punch; the positions of the rigid die wall and the lower punch are fixed; the change in density of the compact is determined by the change in displacement of the upper punch. The downward motion of the upper punch results in a downward force on the particles in contact with it, which in turn impart this force to their neighboring particles in the same direction. This process continues until the forces reach the die wall and the lower punch. It should be noted that the simulation results are not affected by the strain rate since materials are assumed to be strain rate insensitive. Each particle in the composite powder is independent and deformable; the relationship between the positions of particles can be determined in a contact table to judge their contact state in model setup. Upon compaction, the contact between particles was controlled by ‘node to segment method’ which is widely used in FEM simulation. The materials are set to be elastoplastic. Materials parameters for soft and hard particles used in the simulation are listed in Table1. The friction coefficient between particles is set to be 0.2 according to the Coulomb friction model being used. The punches and the die are rigid with smooth surface to reduce their friction with the powder. Once all the simulation parameters are determined, the simulation is run on the server.Table 1
Materials parameters of soft and hard particles used in the simulation.
MaterialsYoung’s modulus Y/GPaPoisson’s ratioν/-Hardening index n/-Strength coefficient C/MPaSoft62.590.330.25225.9Hard1200.250.237575.3


## 2.1. Model Setup

The particle scale MPFEM numerical simulation was carried out using commercial MSC Marc software. In order to reduce constraints and obtain a better simulation of the physical process, DEM dynamic modelling was employed to generate initial random packing structure which was then imported into the FEM model for compaction. The DEM method has been successfully applied to packing of various spherical and nonspherical particulate systems [24–27]. Figure 1(a) shows the DEM generation of initial random packing structure with discretized mesh division on all particles before compaction when the mass fraction of soft particles is Xsoft = 25wt.%, where each particle includes 1014 units and 872 nodes as shown in Figure 1(b). In this case, both soft and hard particles have equal size (1 mm in diameter). At the interface between two particles in contact, the Coulomb friction model was used:(1)σfr=−μσn2πarctanvrrvcnstt,where sfr is the tangential (friction) stress; sn is the normal stress; µ is the friction coefficient; t is the tangential vector in the direction of the relative velocity; vr represents the relative sliding. Physically, the value of rvcnst can be regarded as the value of the relative velocity below which sticking occurs.Figure 1
DEM generated 3D initial random dense packing structure (a) and mesh division of an individual particle in the soft and hard particles (b), the packing structures for with vibration (c) and (d) without vibration.
(a)(b)(c)(d)In the simulation, the yield stress can be expressed by(2)σy=Cε0+ε¯n+Dε¯˙k,where e0 is initial yield strain; ε¯ is equivalent strain; ε¯˙ is equivalent strain rate; parameters C, D, n, and k are material constants. Initially, the equivalent strain and equivalent strain rate were assumed to be zero so that the initial yield strain can be calculated by material constants C and n and Young’s modulus Y in von Mises material model. The von Mises stress is given by(3)σ¯=σ1−σ22+σ2−σ32+σ3−σ121/22,where s1, s2, and s3 are the principal Cauchy stresses along three main axes. The equivalent strain rate is(4)ε¯˙=23εijεij1/2,where eij is the strain tensor.MSC Marc uses Cauchy stress (true stress) and logarithmic strain with updated Lagrange formulation. It is instructive to derive the stiffness matrices for the updated Lagrangian formulation starting from the virtual work principle in the following equation:(5)∫V0SijδEijdV=∫V0bi0δηidV+∫A0ti0δηidA.Here,Sij is the second Piola-Kirchhoff stress; Eij is the Green-Lagrange strain; bi0 is the body force in the reference configuration; ti0 is the traction vector in the reference configuration; and ηi is the virtual displacement. Direct linearization of the left-hand side of (5) yields(6)∫V0SijdδEijdV=∫Vn+1∇ηikσkj∇Δuijdv,where Δu and η are actual incremental and virtual displacements, respectively, and skj is the Cauchy stress tensor. Equation (6) can be further calculated as(7)∫V0dSijδEijdV=∫Vn+1∇sηijLijk1∇sΔuk1dv,where ∇S denotes the symmetric part of ∇, which represents the gradient operator in the present configuration. Lijk1 is convected to the present configuration to yield:(8)K1+K2δu=F−R,where F and R are the external and internal forces, respectively; δu is the correction displacement vector; K1 and K2 are, respectively, the material stiffness matrix and geometric stiffness matrix, given by(9)K1ij=∫Vn+1βimnLmnpqβpqjdv,and(10)K2ij=∫Vn+1σk1Ni,kNj,1dv,where ßimn is the symmetric gradient operator evaluated in the present configuration; Lmnpq represents the material moduli tensor; sk1 is the Cauchy stresses; Ni,k is the shape function gradient matrix.Keeping in view that the reference state is the present state, a rate formulation can be expressed as(11)∫V0S˙ijδEij+Sij∂vk∂Xi∂δηk∂XjdV=∫V0b˙iδηidV+∫A0t˙iδηidA,where S˙ij is the rate of the second Piola-Kirchhoff stress. Equation (11) can be obtained by setting(12)Fij=δij,δEij=δdij,∂∂Xi=∂∂xi,Sij=σij,where Fij is the deformation tensor; dij is the rate of deformation; Xi is the coordinate of initial node; and xi is the coordinate of the node. Hence:(13)∫Vn+1σ∇ijδdij+σij∂vk∂xi∂δηk∂xjdv=∫Vn+1b˙iδηidv+∫Vn+1t˙iδηida,in which b˙i and t˙i are the body force and surface traction in the present configuration, respectively. And in this equation, σij∇ is the Truesdell rate of Cauchy stress.MSC Marc uses Cauchy stress (true stress) and logarithmic strain with updated Lagrange formulation (explained in more detail in the preliminary work [44, 45]).


## 2.2. Simulation Conditions

The boundary conditions of the MPFEM simulation are given as follows: displacement boundary condition is imposed via rigid upper punch; the positions of the rigid die wall and the lower punch are fixed; the change in density of the compact is determined by the change in displacement of the upper punch. The downward motion of the upper punch results in a downward force on the particles in contact with it, which in turn impart this force to their neighboring particles in the same direction. This process continues until the forces reach the die wall and the lower punch. It should be noted that the simulation results are not affected by the strain rate since materials are assumed to be strain rate insensitive. Each particle in the composite powder is independent and deformable; the relationship between the positions of particles can be determined in a contact table to judge their contact state in model setup. Upon compaction, the contact between particles was controlled by ‘node to segment method’ which is widely used in FEM simulation. The materials are set to be elastoplastic. Materials parameters for soft and hard particles used in the simulation are listed in Table1. The friction coefficient between particles is set to be 0.2 according to the Coulomb friction model being used. The punches and the die are rigid with smooth surface to reduce their friction with the powder. Once all the simulation parameters are determined, the simulation is run on the server.Table 1
Materials parameters of soft and hard particles used in the simulation.
MaterialsYoung’s modulus Y/GPaPoisson’s ratioν/-Hardening index n/-Strength coefficient C/MPaSoft62.590.330.25225.9Hard1200.250.237575.3


## 3. Results and Discussion

### 3.1. Compaction Process and Macroscopic Analyses

In the forming process, the influence of the initial packing structure of the powder was considered in its entirety. In addition to the initial random packing before compaction, in order to increase the loose density of the composite powder, mechanical vibration to the composite powder in the DEM simulation was applied. The forming process of these powders was studied by numerical simulation and physical experiment, as shown in Figure2. It can be seen that the results of physical experiment and numerical simulation are in agreement, which proves the effectiveness of the model used in our numerical simulation. In addition, we can also find that the change law of the three-dimensional compaction curve of the composite powder is similar to our previous two-dimensional simulation results [40, 44], and the specific reasons will not be described here. It should be mentioned that, during the formation of both soft and hard particles, although the initial packing structure of powder has little effect on the relative density of final compacts, it can be seen from Figure 2 that it has a significant influence on the early stage of forming process. At the same time, it also has a significant impact on the microproperties of the compaction, which will be discussed in the following sections.Figure 2
Numerical and physical compaction curves for the compaction of initial packing structures with and without mechanical vibration.In contrast to the macroscale continuous simulation method, our MPFEM simulation can also characterize the morphology and internal stress of the compact during the whole compaction process from the particle scale, as shown in Figure3. In this case, the soft particles content in the composite powder is 25wt%; vibration treatment is applied to the powder before compaction. It can be seen from the figure that, in contrast to the two-dimensional compaction process, the evolution process of yield surface can be clearly observed in the three-dimensional compaction process of composite particles. This means that the mechanism of powder compaction densification can be accurately simulated from the microscopic point of view. As the compaction process evolves, it is found that the deformations of both soft and hard particles in the compaction gradually increase. In particular, the deformation of the soft particles is greater than that of the hard particles, and the stress concentration is mainly generated in the hard particles, which is similar to the stress concentration characteristics in two-dimensional compaction process [36]. It should be noted that the deformation of soft and hard particles gradually changes from curved surface to plane at the part contacting with the die wall, where the deformation amount of the two particles is almost equal, although the stress concentration of the deformation part of hard particles is still high. However, due to the release of internal stress after the pressure is relieved, the residual stress of both soft and hard particles near the edge of the die wall is almost equal, except for the residual stress in the hard particles that is evenly distributed. In addition, the size is almost equal; however, only the soft particles are under higher residual stress in the deformed part. This phenomenon is difficult to reflect in the two-dimensional forming process of composite powder [46, 47].Figure 3
Evolution of the morphology and internal equivalent von Mises stress distribution of composite compact from loading to pressure releasing with mechanical vibration when (X)soft = 25wt%. (a)–(f) correspond to (P) = 7, 50, 100, and 200 MPa, pressure holding, and pressure releasing, respectively.In order to further explore the mechanism that governs the filling of internal pores of the compacts and the evolution of the particle morphology during the three-dimensional compaction of composite powders, the compaction obtained under the condition of pressureP = 200 MPa was layered and displayed, and the compaction morphology of different heights and the corresponding stress distribution were obtained as shown in Figure 4. It can be seen from the figure that the density at the bottom of the compact is lower than that at the top, which is caused by the top-down unidirectional die pressing used in this study. Since the force is transmitted from top to bottom through the upper die punch, the powder near the upper die punch is subjected to greater pressure and the degree of deformation during the forming process is relatively high. However, this deformation is not uniform. Usually, the pores in the contact parts of soft particles are filled more than the pores in the contact parts of hard particles. This conclusion is in good agreement with the maximum and minimum density distribution at the top edge and bottom edge of the compact obtained by continuous macroscopic simulation [20, 22]. In addition, it can be seen from the figure that the deformation of soft particles is far greater than that of hard particles, and the closer the soft particles are to the upper die, the greater the deformation of soft particles. In addition, the surface is flatter, and the deformation of the lower soft particles is mainly caused by the action of hard particles. The surface of the soft particle that is in contact with the hard particle is inwardly curved. The irregular particle shape produced by this depression is difficult to squeeze into the adjacent pores inside the compact, which hinders further filling of the pores.Figure 4
Equivalent von Mises stress distributions for local structures obtained from different heights of the compact when (P) = 200 MPa. (a)–(c) correspond to upper part, middle part, and bottom part of the compact.


### 3.2. Microscopic Property Characterization and Densification Mechanism Analysis

In addition to the macroscopic properties, we also characterized and compared the microscopic properties of the composite powder during the forming process from the particle scale. Figure5 shows the morphology and stress changes of single hard particle and contact area of soft particle in the center of the compaction at different stages of the pressing process of both soft and hard particles, where different colors of particles represent stress levels. Three states in the process of compaction were selected for study, and the relative densities of the compaction were 0.64, 0.74, and 0.84, respectively. As can be seen from the figure, with the increase of pressing pressure, hard particles barely deform, but their color changes are manifest, which indicates that hard particles contact with other particles more and more closely, and the interaction between particles is also increasingly intense. On the other hand, the soft particles underwent large deformations, and the parts with significant deformation even exhibit inward depression, while the parts with bulges show almost no color change. This indicates that the part has not been in contact with other particles during the forming process or that the force in the contact area is very small even when it is in contact. The filling of pores is accomplished by the mass transfer at the protrusion site of soft particles. The stress at the deformation site where soft particles are in close contact with other particles is almost equal. The specific changes and their causes will be discussed systematically in the next part of the quantitative analysis. In addition, the 3D simulation results obtained here are similar to the 2D simulation results [44].Figure 5
Evolution of the morphology of individual particles at different stages of compaction of vibrated initial packing structure of composite powder.In order to further quantitatively characterize the deformation and mechanical behavior of the powder particles in the forming process, both soft and hard particles presented in Figure4 were studied in depth, focusing on the distribution of internal stress in the two different particles and the evolution of residual stress after unloading as shown in Figure 6. As can be seen from the figure, the stress value in the center of hard particles in the compacted blank formed at P = 200 MPa is far less than that at the surface with a large deformation, which is consistent with the two-dimensional simulation results. As can be seen from Figure 6, the locations with comparatively large particle shape variable correspond to locations where other hard particles contact each other. Since the yield strength of hard particles is higher than that of soft particles, the plastic deformation of hard particles requires more external energy, which causes greater concentration of stress at the plastic deformation site. However, for the soft particles in contact with it, the stress generated in the central part is not significantly different from the part with larger plastic deformation. There are two kinds of local deformation morphology of soft particles; one tends to be flat, and the other is concave to the inside of soft particles. The morphology that tends to be flat indicates that it is situated at the contact of two soft particles. As the pressing progresses, the contact between the two particles changes from point to surface. The inwardly concave morphology of soft particles indicates that the contact is between soft and hard particles. Due to the high hardness of hard particles, the deformation of hard particles is very small during the forming process. Under the action of the hard particles, soft particles exhibit large plastic deformation and can envelop the hard particles. In this case, the maximum stress is concentrated in the depression of soft particles. After the pressure is relieved, the internal stress of the particles is redistributed and the self-equilibrium of stress is achieved. It should be noted that the stress in the center of both soft and hard particles is reduced almost to zero, and comparatively high residual stress remains only at the edges and corners of the particles.Figure 6
Equivalent von Mises stress distributions and morphology evolution of particles in Figure5 at (P) = 200 MPa and unloading. (a) and (b) show a hard particle and (c) and (d) show a soft particle; left figures correspond to (P) = 200 MPa and right to unloading.The stress distribution inside hard particles in Figure6 was quantitatively characterized. Figure 7(a) shows the equivalent von Mises stress distribution inside hard particles along three axes from the center (0.5 mm). It can be seen from the figure that the equivalent von Mises stress value of the hard particle varies from 175 MPa to 340 MPa and the distribution is not uniform. The stress value in the center of the particle is about 230 MPa, and it fluctuates greatly in the process of extending to the edge of the particle. However, regardless of the axis of direction, greater equivalent von Mises stress values are found closer to the edge of the particle. In the process of unidirectional die pressing, the upper die moves along the negative direction of Z axis, and the maximum stress on the upper surface of hard particles is about 340 MPa. Figure 7(b) shows the distribution of equivalent von Mises stress in the center of hard particle (0.5 mm) along three axes after pressure relief. It can be seen that the equivalent von Mises stress in the center of the hard particles drops from 230 MPa to 80 MPa after the pressure is relieved, and the distribution from the center to both sides is more uneven. The equivalent von Mises stress value varied from 60 MPa to 250 MPa. Although it decreased after pressure relief, the difference between maximum and minimum is very close the equivalent difference inside the particle at P = 200 MPa. The larger internal stress difference is not conducive to the subsequent sintering process.Figure 7
Equivalent von Mises stress distributions of hard particle (from center (0.5 mm) to outside) along (X)Y and Z directions. (a) P = 200 MPa; (b) unloading.
(a)(b)In a similar fashion, we also quantitatively characterized the internal stress distribution of the soft particles in Figure6. Figure 8(a) shows the equivalent von Mises stress distribution from the center (0.5 mm) of the soft particles along the three coordinate axes. It can be seen from the figure that the equivalent von Mises stress value of the soft particles varies between 137 MPa and 240 MPa, and the internal stress distribution is more uneven than that of the hard particles in Figure 7(a). This is because the deformation of the soft particles is greater than that of the hard particles during the forming process, which generates multiple shear forces inside the soft particles. The stress value in the center of soft particles is about 175 MPa, which is 55 MPa lower than that in the center of hard particles. Although the upper die moves along the negative direction of the Z axis, the equivalent von Mises stress at the lower part of the particle is the highest. This is because the lower part of the particle is influenced by the action of hard particles (the most severe part of the depression in Figure 5 is the lower part of the soft particle), and the maximum stress on the lower surface of the particle is about 240 MPa. Figure 8(b) shows the distribution of equivalent von Mises stress along three axes from the soft particle center (0.5 mm) after pressure relief. By comparing the equivalent von Mises stress distribution within hard particles after pressure relief (Figure 7(b)), it can be seen that the stress distribution within soft particles is more uneven. However, the equivalent von Mises stress value within soft particles varies between 25 MPa and 175 MPa, and the difference is smaller than that within hard particles after pressure relief. Therefore, in order to obtain high quality compaction, we should devote ourselves to relieving the larger stress gradient inside the particle and reducing the stress gradient inside the particle.Figure 8
Equivalent von Mises stress distributions in soft particle from center (0.5 mm) to exterior along (X)Y and Z directions. (a) (P) = 200 MPa, (b) unloading.
(a)(b)In order to further study the characteristics of the maximum stress generated on the surface of both soft and hard particles, the ones in the center of the compaction were considered as the objects, and the points producing the maximum stress on their surfaces were selected, among which the node on soft particle was Node17672 and the node on hard particle was Node52206. In addition, both soft and hard particles near the die wall were also considered as objects, and the points producing the maximum stress value on their surfaces were selected, among which there is the node Node62086 on soft particle and the node Node23250 on hard particle. The stress present at these nodes was compared and analyzed. It should be noted that these points represent the maximum stress points in the compaction, and the selected particles are the center of compaction and close to the die wall, in order to be representative and universal. The results are shown in Figure9. As can be seen from the figure, the equivalent von Mises stress values at these nodes increase continuously with the process of formation. However, the equivalent von Mises stress values of these nodes decrease after pressure relief. By contrast, the equivalent von Mises stress value at soft particle nodes is much smaller than that at hard particle nodes, and the difference becomes larger with the increase of the number of simulated steps. The maximum difference of equivalent von Mises stress between soft and hard particles in the center at the same number of simulation steps is about 200 MPa, and the maximum difference near the die wall with the same number of simulation steps is about 300 MPa. However, the equivalent von Mises stress values for the soft particles in the center and near the die wall are not significantly different, while the equivalent von Mises stress values of the hard particles in the center and near the die wall are significantly different with the maximum difference of 220 MPa. This is because the deformation of soft particles in any position is very large during the forming process, while the hard particles near the die wall exhibit greater deformation than in the center. After the pressure is relieved, the equivalent von Mises stress value of the particle near the die wall decreased more than that of the particle at the center. However, the residual stress at the particles near the die wall remains comparatively large, node17672 about 187 MPa; node52206 about 337 MPa; node62086 about 80 MPa; node23250 about 385 MPa.Figure 9
Variation of equivalent von Mises stress with incremental modelling step for different particle nodes in the compact during compaction.


## 3.1. Compaction Process and Macroscopic Analyses

In the forming process, the influence of the initial packing structure of the powder was considered in its entirety. In addition to the initial random packing before compaction, in order to increase the loose density of the composite powder, mechanical vibration to the composite powder in the DEM simulation was applied. The forming process of these powders was studied by numerical simulation and physical experiment, as shown in Figure2. It can be seen that the results of physical experiment and numerical simulation are in agreement, which proves the effectiveness of the model used in our numerical simulation. In addition, we can also find that the change law of the three-dimensional compaction curve of the composite powder is similar to our previous two-dimensional simulation results [40, 44], and the specific reasons will not be described here. It should be mentioned that, during the formation of both soft and hard particles, although the initial packing structure of powder has little effect on the relative density of final compacts, it can be seen from Figure 2 that it has a significant influence on the early stage of forming process. At the same time, it also has a significant impact on the microproperties of the compaction, which will be discussed in the following sections.Figure 2
Numerical and physical compaction curves for the compaction of initial packing structures with and without mechanical vibration.In contrast to the macroscale continuous simulation method, our MPFEM simulation can also characterize the morphology and internal stress of the compact during the whole compaction process from the particle scale, as shown in Figure3. In this case, the soft particles content in the composite powder is 25wt%; vibration treatment is applied to the powder before compaction. It can be seen from the figure that, in contrast to the two-dimensional compaction process, the evolution process of yield surface can be clearly observed in the three-dimensional compaction process of composite particles. This means that the mechanism of powder compaction densification can be accurately simulated from the microscopic point of view. As the compaction process evolves, it is found that the deformations of both soft and hard particles in the compaction gradually increase. In particular, the deformation of the soft particles is greater than that of the hard particles, and the stress concentration is mainly generated in the hard particles, which is similar to the stress concentration characteristics in two-dimensional compaction process [36]. It should be noted that the deformation of soft and hard particles gradually changes from curved surface to plane at the part contacting with the die wall, where the deformation amount of the two particles is almost equal, although the stress concentration of the deformation part of hard particles is still high. However, due to the release of internal stress after the pressure is relieved, the residual stress of both soft and hard particles near the edge of the die wall is almost equal, except for the residual stress in the hard particles that is evenly distributed. In addition, the size is almost equal; however, only the soft particles are under higher residual stress in the deformed part. This phenomenon is difficult to reflect in the two-dimensional forming process of composite powder [46, 47].Figure 3
Evolution of the morphology and internal equivalent von Mises stress distribution of composite compact from loading to pressure releasing with mechanical vibration when (X)soft = 25wt%. (a)–(f) correspond to (P) = 7, 50, 100, and 200 MPa, pressure holding, and pressure releasing, respectively.In order to further explore the mechanism that governs the filling of internal pores of the compacts and the evolution of the particle morphology during the three-dimensional compaction of composite powders, the compaction obtained under the condition of pressureP = 200 MPa was layered and displayed, and the compaction morphology of different heights and the corresponding stress distribution were obtained as shown in Figure 4. It can be seen from the figure that the density at the bottom of the compact is lower than that at the top, which is caused by the top-down unidirectional die pressing used in this study. Since the force is transmitted from top to bottom through the upper die punch, the powder near the upper die punch is subjected to greater pressure and the degree of deformation during the forming process is relatively high. However, this deformation is not uniform. Usually, the pores in the contact parts of soft particles are filled more than the pores in the contact parts of hard particles. This conclusion is in good agreement with the maximum and minimum density distribution at the top edge and bottom edge of the compact obtained by continuous macroscopic simulation [20, 22]. In addition, it can be seen from the figure that the deformation of soft particles is far greater than that of hard particles, and the closer the soft particles are to the upper die, the greater the deformation of soft particles. In addition, the surface is flatter, and the deformation of the lower soft particles is mainly caused by the action of hard particles. The surface of the soft particle that is in contact with the hard particle is inwardly curved. The irregular particle shape produced by this depression is difficult to squeeze into the adjacent pores inside the compact, which hinders further filling of the pores.Figure 4
Equivalent von Mises stress distributions for local structures obtained from different heights of the compact when (P) = 200 MPa. (a)–(c) correspond to upper part, middle part, and bottom part of the compact.


## 3.2. Microscopic Property Characterization and Densification Mechanism Analysis

In addition to the macroscopic properties, we also characterized and compared the microscopic properties of the composite powder during the forming process from the particle scale. Figure5 shows the morphology and stress changes of single hard particle and contact area of soft particle in the center of the compaction at different stages of the pressing process of both soft and hard particles, where different colors of particles represent stress levels. Three states in the process of compaction were selected for study, and the relative densities of the compaction were 0.64, 0.74, and 0.84, respectively. As can be seen from the figure, with the increase of pressing pressure, hard particles barely deform, but their color changes are manifest, which indicates that hard particles contact with other particles more and more closely, and the interaction between particles is also increasingly intense. On the other hand, the soft particles underwent large deformations, and the parts with significant deformation even exhibit inward depression, while the parts with bulges show almost no color change. This indicates that the part has not been in contact with other particles during the forming process or that the force in the contact area is very small even when it is in contact. The filling of pores is accomplished by the mass transfer at the protrusion site of soft particles. The stress at the deformation site where soft particles are in close contact with other particles is almost equal. The specific changes and their causes will be discussed systematically in the next part of the quantitative analysis. In addition, the 3D simulation results obtained here are similar to the 2D simulation results [44].Figure 5
Evolution of the morphology of individual particles at different stages of compaction of vibrated initial packing structure of composite powder.In order to further quantitatively characterize the deformation and mechanical behavior of the powder particles in the forming process, both soft and hard particles presented in Figure4 were studied in depth, focusing on the distribution of internal stress in the two different particles and the evolution of residual stress after unloading as shown in Figure 6. As can be seen from the figure, the stress value in the center of hard particles in the compacted blank formed at P = 200 MPa is far less than that at the surface with a large deformation, which is consistent with the two-dimensional simulation results. As can be seen from Figure 6, the locations with comparatively large particle shape variable correspond to locations where other hard particles contact each other. Since the yield strength of hard particles is higher than that of soft particles, the plastic deformation of hard particles requires more external energy, which causes greater concentration of stress at the plastic deformation site. However, for the soft particles in contact with it, the stress generated in the central part is not significantly different from the part with larger plastic deformation. There are two kinds of local deformation morphology of soft particles; one tends to be flat, and the other is concave to the inside of soft particles. The morphology that tends to be flat indicates that it is situated at the contact of two soft particles. As the pressing progresses, the contact between the two particles changes from point to surface. The inwardly concave morphology of soft particles indicates that the contact is between soft and hard particles. Due to the high hardness of hard particles, the deformation of hard particles is very small during the forming process. Under the action of the hard particles, soft particles exhibit large plastic deformation and can envelop the hard particles. In this case, the maximum stress is concentrated in the depression of soft particles. After the pressure is relieved, the internal stress of the particles is redistributed and the self-equilibrium of stress is achieved. It should be noted that the stress in the center of both soft and hard particles is reduced almost to zero, and comparatively high residual stress remains only at the edges and corners of the particles.Figure 6
Equivalent von Mises stress distributions and morphology evolution of particles in Figure5 at (P) = 200 MPa and unloading. (a) and (b) show a hard particle and (c) and (d) show a soft particle; left figures correspond to (P) = 200 MPa and right to unloading.The stress distribution inside hard particles in Figure6 was quantitatively characterized. Figure 7(a) shows the equivalent von Mises stress distribution inside hard particles along three axes from the center (0.5 mm). It can be seen from the figure that the equivalent von Mises stress value of the hard particle varies from 175 MPa to 340 MPa and the distribution is not uniform. The stress value in the center of the particle is about 230 MPa, and it fluctuates greatly in the process of extending to the edge of the particle. However, regardless of the axis of direction, greater equivalent von Mises stress values are found closer to the edge of the particle. In the process of unidirectional die pressing, the upper die moves along the negative direction of Z axis, and the maximum stress on the upper surface of hard particles is about 340 MPa. Figure 7(b) shows the distribution of equivalent von Mises stress in the center of hard particle (0.5 mm) along three axes after pressure relief. It can be seen that the equivalent von Mises stress in the center of the hard particles drops from 230 MPa to 80 MPa after the pressure is relieved, and the distribution from the center to both sides is more uneven. The equivalent von Mises stress value varied from 60 MPa to 250 MPa. Although it decreased after pressure relief, the difference between maximum and minimum is very close the equivalent difference inside the particle at P = 200 MPa. The larger internal stress difference is not conducive to the subsequent sintering process.Figure 7
Equivalent von Mises stress distributions of hard particle (from center (0.5 mm) to outside) along (X)Y and Z directions. (a) P = 200 MPa; (b) unloading.
(a)(b)In a similar fashion, we also quantitatively characterized the internal stress distribution of the soft particles in Figure6. Figure 8(a) shows the equivalent von Mises stress distribution from the center (0.5 mm) of the soft particles along the three coordinate axes. It can be seen from the figure that the equivalent von Mises stress value of the soft particles varies between 137 MPa and 240 MPa, and the internal stress distribution is more uneven than that of the hard particles in Figure 7(a). This is because the deformation of the soft particles is greater than that of the hard particles during the forming process, which generates multiple shear forces inside the soft particles. The stress value in the center of soft particles is about 175 MPa, which is 55 MPa lower than that in the center of hard particles. Although the upper die moves along the negative direction of the Z axis, the equivalent von Mises stress at the lower part of the particle is the highest. This is because the lower part of the particle is influenced by the action of hard particles (the most severe part of the depression in Figure 5 is the lower part of the soft particle), and the maximum stress on the lower surface of the particle is about 240 MPa. Figure 8(b) shows the distribution of equivalent von Mises stress along three axes from the soft particle center (0.5 mm) after pressure relief. By comparing the equivalent von Mises stress distribution within hard particles after pressure relief (Figure 7(b)), it can be seen that the stress distribution within soft particles is more uneven. However, the equivalent von Mises stress value within soft particles varies between 25 MPa and 175 MPa, and the difference is smaller than that within hard particles after pressure relief. Therefore, in order to obtain high quality compaction, we should devote ourselves to relieving the larger stress gradient inside the particle and reducing the stress gradient inside the particle.Figure 8
Equivalent von Mises stress distributions in soft particle from center (0.5 mm) to exterior along (X)Y and Z directions. (a) (P) = 200 MPa, (b) unloading.
(a)(b)In order to further study the characteristics of the maximum stress generated on the surface of both soft and hard particles, the ones in the center of the compaction were considered as the objects, and the points producing the maximum stress on their surfaces were selected, among which the node on soft particle was Node17672 and the node on hard particle was Node52206. In addition, both soft and hard particles near the die wall were also considered as objects, and the points producing the maximum stress value on their surfaces were selected, among which there is the node Node62086 on soft particle and the node Node23250 on hard particle. The stress present at these nodes was compared and analyzed. It should be noted that these points represent the maximum stress points in the compaction, and the selected particles are the center of compaction and close to the die wall, in order to be representative and universal. The results are shown in Figure9. As can be seen from the figure, the equivalent von Mises stress values at these nodes increase continuously with the process of formation. However, the equivalent von Mises stress values of these nodes decrease after pressure relief. By contrast, the equivalent von Mises stress value at soft particle nodes is much smaller than that at hard particle nodes, and the difference becomes larger with the increase of the number of simulated steps. The maximum difference of equivalent von Mises stress between soft and hard particles in the center at the same number of simulation steps is about 200 MPa, and the maximum difference near the die wall with the same number of simulation steps is about 300 MPa. However, the equivalent von Mises stress values for the soft particles in the center and near the die wall are not significantly different, while the equivalent von Mises stress values of the hard particles in the center and near the die wall are significantly different with the maximum difference of 220 MPa. This is because the deformation of soft particles in any position is very large during the forming process, while the hard particles near the die wall exhibit greater deformation than in the center. After the pressure is relieved, the equivalent von Mises stress value of the particle near the die wall decreased more than that of the particle at the center. However, the residual stress at the particles near the die wall remains comparatively large, node17672 about 187 MPa; node52206 about 337 MPa; node62086 about 80 MPa; node23250 about 385 MPa.Figure 9
Variation of equivalent von Mises stress with incremental modelling step for different particle nodes in the compact during compaction.


## 4. Conclusions

In this paper, the MPFEM method was used to simulate the three-dimensional forming process of composite powder at the particle scale. The effect of the initial packing structure of the powder before compaction was investigated. The macroscopic and microscopic properties of the composite powders during the compaction process were quantitatively characterized and analyzed. In addition, the densification mechanism of the composite powders was studied from the perspective of pore filling and force. The main conclusions obtained are as follows:(1)
In the 3D numerical simulation of the forming process of composite powders, the process of pressing after vibration is applied to the powder before forming can result in higher relative density under low pressure. The vibration can result in more favorable arrangement of the particles and eliminate the bridge and arch structures that hinder the uniform distribution of pores in the initial packing structure.(2)
The 3D numerical simulation can reproduce the pore filling mechanism of the composite powder compaction, evolution of the internal stress, and particle morphology of the composite powder compact in real time. After vibration is applied to the powder before forming, the distribution of pores in the compacted blank after forming is more uniform than that without vibration, and the pores close to the upper die part of the compacted blank are smaller than those in the lower part of the compacted compact. The surface deformation of the particles near the die wall and the upper and lower parts of the punching die is relatively smooth, while the interior particles, particularly the soft particles, exhibit concave surface due to the extrusion of hard particles. The maximum stress occurs at the edges and corners of deformation and in the center of plane and concave with severe deformation.(3)
For either soft or hard particles, the stress is concentrated at their central part during the process of pressing, and the internal stress distribution is not uniform. After the pressure is relieved, the stress concentration in the center of the particle is greatly reduced, and the stress is rebalanced. The larger stress concentration appears at the edges and corners of the deformed particle, which causes the internal stress distribution of the particle after the pressure is relieved to be extremely uneven. In contrast, the stress inhomogeneity inside soft particles is higher than that inside hard particles. In the subsequent sintering process, cracks occur not only at the contact points of the particles, but also at the edges and corners of the particles, resulting in a decrease in the quality of sintered compact.(4)
The results of the physical experiment are in good agreement with the numerical simulation, which proves the validity of the numerical model used in the MPFEM simulation. The results of physical experiments show that applying vibration to the powder before forming can increase the loose packing density of the composite powder, and simultaneously time improves both macro- and microproperties of the compact.



---

*Source: 8580849-2022-02-10.xml*