Influence of particle deformation on the plastic flow of ductile granular materials

A multi-particle finite-element method was proposed to study the elastic-plastic behaviour of ductile powders composed of highly deformable elastic-plastic particles. The focus was put on the study of the uniqueness of the direction of plastic strain increment vectors for a given stress state on the plastic limit, which was assessed using a spherical stress-probing method. Results revealed a non-uniqueness of the direction of plastic flow in a small region of the stress space located in the vicinity of the loading point. The direction of plastic flow was almost unique elsewhere on the plastic limit. The non-uniqueness was explained using a combination of two distinct mechanisms for plastic deformation involving two very different plastic limits.


Introduction
This paper is concerned with the mechanical behaviour of granular assemblies involving highly deformable elasticplastic particles. Such assemblies will be referred to as 'ductile powders'. The main targeted application is powder compaction. There is a need for continuummechanics-based constitutive models which would be able to predict the evolution of stress and strain with a reasonable accuracy in localised high-shear zones where cracking is susceptible to happen.
To this end, a numerical model involving an assembly of meshed spheres interacting through contact conditions is studied in the framework of the finite-element method; this approach will be named Multi-Particle Finite-Element Method (MPFEM) in the following. This method allows for the simulation of many loading paths that are not attainable to experimental studies due to technical difficulties. By studying the response of the numerical sample (which is considered as a model material for ductile powders) to such loading paths, insight is given on the mechanical behaviour of such materials.
Previous studies [1][2][3] showed the ability of the MPFEM in deriving yield surfaces and describing strainhardening mechanisms and their relation to the evolution of the microstructure. The focus of the present work is the study of plastic flow.

Numerical model
In this study, a polydisperse assembly of 50 randomlydistributed spherical particles was built into the commercial finite-element software ABAQUS using an explicit e-mail: barthelemy.harthong@3sr-grenoble.fr time integration scheme. Let the dimensions of the box be taken as 1x1x1 mm, then the stresses presented in the following figures are expressed in MPa. Particles were meshed using approximately 3300 quadratic tetrahedral elements per particle, which allowed to keep a reasonable accuracy in the description of contact surfaces. The effect of the number of particles was studied and results showed that yield surfaces were decreasing in size when increasing the number of particles but were keeping the same shape and orientation. The influence of polydispersity was not studied. It is important to understand that the use of only 50 spheres in the numerical model implies that the following results do not mean to be quantitative, general results for assemblies of ductile particles. The purpose of the present paper is solely to observe the behaviour of the 50spheres sample and provide explanations of this behaviour. Only the understanding of the observed behaviour can be transferred to real materials.
Individual particles were assigned a constitutive model associating linear, isotropic elasticity (with elastic modulus E = 10 GPa and Poisson coefficient ν = 0.435) and Mises plasticity with an isotropic power-law strainhardening: where k = 15.5 MPa and n = 0.35 are material parameters, and σ 0 , σ Y and ε pl are, respectively, the initial and current (Mises) yield stress and the equivalent plastic strain. These parameters were calibrated on a lead alloy corresponding to the experimental tests performed by Chen [4]. Particle/particle and particle/plane contact interactions were modelled by means of a penalty contact algorithm and a classical Coulomb friction law with a friction coefficient f = 0.1. Chen [4] observed that the forces on the planes increased by approximately 10% when raising the particle/plane friction coefficient from 0 to 0.1 at a relative density D = 0.97 (much less than 10% for lower values of D); and another 10% when raising it from 0.1 to 0.25. No contact cohesion was introduced. The particles were enclosed in a cubic box made of six rigid planes which were allowed to move along their normal, so that initially parallel planes remained parallel and initially perpendicular planes remained perpendicular throughout the loading. As a result, distortional strains were all kept to zero. Using such boundary conditions, it was checked that tangent force components on the rigid planes were negligible compared to normal components, meaning that no rotation of either stress or strain tensors was applied. As a result, the numerical sample being initially isotropic, stress and strain principal directions were assumed to remain coaxial throughout the simulation, so that the present analysis is restricted to the 3D space of principal stresses or strains, in the particular case of coaxiality of the principal axes of stress and strain.
The numerical assembly was compacted using either isotropic straining (referred to as isotropic compaction) and uniaxial straining along x direction (referred to as closed-die compaction). Boundary conditions were either stress-driven (the total normal force on each plane was controlled to reach desired stress values) or straindriven (normal displacements were imposed to reach desired strains).

Probing of yield surfaces
The numerical sample was loaded following the isotropic or closed-die loading path and unloaded. From the unloaded state, a series of reloading steps were applied in every direction of the stress space to probe yield surfaces ( Fig. 1).
Following [2], the detection of yield points during the reloading steps was based on the dissipated energy W, equal to the sum of plastic and frictional dissipation within the numerical model. The onset of plasticity was defined as the point where the total dissipation in the assembly reached a threshold value equal to W 0 +0.3%W 0 , W 0 being the value of the total dissipation at the end of the unloading step.   2 shows the yield surfaces in a deviatoric stress/mean stress diagram with various values of final relative density D for the 50-particles numerical sample which was submitted to isotropic and closed die compaction.

Spherical stress probing method
In classical elasto-plasticity, the incremental plastic strain obeys the flow rule equation: where dε pl is the increment of plastic strain, σ the stress tensor, χ the strain-hardening variables, g the plastic potential and λ the plastic multiplier. Assuming g is regular, equation (2) states that for a given state of stress and strain-hardening, the direction of the incremental plastic strain dε pl is unique.
In the present work, the existence of a plastic potential, and more generally, the incremental plastic flow could be studied on the numerical sample by using a spherical stress probing method [5]. This method consisted firstly in loading the numerical sample to a given state of stress σ 0 on the yield surface (i.e. at the elastic limit). Secondly, small stress increments dσ are applied in every direction of the stress space (Fig. 3a). Thirdly, the sample is unloaded back to σ 0 . Assuming that no plastic dissipation occurs during unloading, the incremental plastic strain dε pl was assimilated to the difference in total strain after loading and unloading.
The plastic strain response envelope is the curve formed by the tips of the incremental plastic strain vectors corresponding to the response of the sample to the probing stress increments dσ (Fig. 3b). In the present paper, only stress increments in the Rendulic plane were applied. Referring back to equation (2), if the flow rule applies and if a regular potential function exists, the direction of the plastic strain increment dε pl does not depend on the direction of the stress increment dσ. In such a case, the plastic strain response envelope is a line segment, since the direction of dε pl is unique and entirely defined by σ and χ. If the response envelope of the plastic strain is not a straight line, then dε pl depends on the direction of dσ, which means that the flow rule equation (2) is not valid.

Results and discussion
To improve readability of the results, plastic strain response envelopes were plotted in deviatoric strain/volumetric strain graphs, then normalised and superimposed to deviatoric stress/mean stress graphs, following the common representation. Fig. 4 shows normalised plastic strain response envelopes for isotropic and closed die compaction of the sample for a constant stress increment of magnitude dσ in the case of a relative density D = 0.7. The value of dσ was chosen as small as possible, but large enough so that the incremental plastic strain response was significantly larger than numerical noise. Fig. 4 show that most plastic strain increment vectors are relatively close to line segments, and are approximately normal to the yield surface, suggesting the validity of an associated flow rule for such materials.
However, in the vicinity of the tip of the yield surfaces, the direction of plastic strain increments was far from being unique. Close to the tips of closed-die-compaction yield surfaces, the plastic strain response envelopes were oriented toward a direction which significantly differed The results indicated a correlation between the nonuniqueness of the direction of plastic strain increments and the particular stress state corresponding to the loading point. As pointed out in [6] and [3], the stress state corresponding to the loading point is actually the stress state for which all contact surfaces between particles are approximately normal to contact forces, because it is the stress state which created the contact surfaces (at least for monotonous loadings such as isotropic or closed-die compaction). In other words, the loading point is the state of stress for which shear stresses at contacts are minimum in average.
As a consequence, sliding between particles was less likely to happen for the stress states close to the loading point than for other stress states which mobilise more shear stresses at contact surfaces. This resulted in the fact that the ratio of frictional dissipation over plastic dissipation was minimum at the loading point, as can be seen in Fig. 5, and indicating that for such a stress state, particles were blocked. As a result, the macroscopic plastic strain was mainly, if not solely due to plastic deformation of particles.
Furthermore, the contact surfaces which increased in size during densification opposed an increasing resistance to normal stresses. This phenomenon results in the fact that the sample develops an increased resistance to the stress state of the loading point compared to other stress states. In the present case, this phenomenon resulted in plastic strain increment vectors being very small in magnitude (which was observed even though not shown here).
A plausible explanation for the non-uniqueness of the direction of plastic flow can be proposed on the basis of a simplified analogy. The analogy consists in a block sliding on a plane under the effect of a force F (Fig. 6). Let F increase in magnitude until the block reaches the limit of sliding (the analogous of the yield point). At this stage, let a small force dF be added in an arbitrary direction (the analogous of the probing stress increment). For some di- rections of dF, the block will pass the sliding limit and move on the plane, resulting in a small displacement du (the analogous of the plastic strain increment vector). In such an analogy, the block represents an average particle, and the plane is the contact structure which defines the possible directions in space along which the particle is allowed to move. The force F represents the stress state. The motion of the block will be the result of 1-the direction of the total force F + dF and 2-the orientation of the plane. dF being very small, the approximation that F + dF = F is valid and it follows that du is uniquely defined by F. In other words, the block moves in a direction defined by the projection of F on the plane (Fig. 6a).
However, the previous reasoning is no more valid when the force F is normal to the plane. In such a case, if dF remains extremely small, it is not possible to reach the sliding limit. However it is possible to imagine that if the friction coefficient is low enough, dF can be large enough to make the block pass the sliding limit while still very small compared to F. In such a case, the direction of slip, du, is purely determined by the direction of dF since the projection of F on the plane is zero. This result implies that the direction of slip is non-unique: it depends on dF and not solely on F.
To improve the analogy, it is possible to introduce a deformable plane such that a very small component of displacement is allowed in the direction of the normal to the plane, but this component is so small that in the case where F is normal to the plane, it remains comparable to the component due to dF alongside the plane.
The sliding block analogy suggests that the requirement to observe the non-uniqueness of the direction of plastic flow is mainly to have two deformation mechanisms with very different plastic limits (in the analogy, the movement against the plane and alongside the plane). In the case of the 50-particle numerical sample studied here, the two mechanisms at stake were thought to be the shearing and normal loading at contact surfaces. As mentioned above, the direction of the stress space corresponding to the loading path is the one for which the sample shows the best resistance. Let this direction be called a "strong" direction, meaning a direction in which contact surfaces are the largest in average, while other stress directions are "weak". When an average particle is loaded in a "strong" direction, the corresponding deformation or displacement of the particle is very small, whereas it can be significantly deformed by a relatively small stress when loaded in a "weak" direction. The "strong" and "weak" directions are equivalent to the directions normal and parallel to the supporting plane in the analogy. In such a case, the direction of plastic strain is influenced by the direction of the stress increment, hence the non-unique behaviour.

Conclusion
A numerical assembly of 50 deformable particles with elastic-plastic behaviour was studied with the aim of characterising incremental plastic flow. Results show that the direction of plastic flow exhibited a well-marked nonuniqueness for some stress states. The non-uniqueness was correlated with the fact that particles displacements were restrained and stresses at contacts were mainly normal. These results suggested the existence of to distinct mechanisms for macroscopic plastic strains in such materials, so that a very high stress in a particular direction could lead to a very small plastic strain while very small stresses in other directions could lead to plastic strains of similar magnitude, leading to the non-uniqueness of plastic flow direction.