Modeling strain and density distributions during high-pressure torsion of pre-compacted powder materials

ABSTRACT It is shown by finite element modeling that under high-pressure torsion of water-atomized pure iron powder a zone with high shear strain and high density develops first in a layer adjacent to the moving anvil. The finite element simulation results and an analytical calculation using the minimum plastic energy extremum principle both show that with further deformation the compacted zone extends progressively in the radial and axial directions throughout the specimen. GRAPHICAL ABSTRACT IMPACT STATEMENT The strain non-uniformities that appear during compaction of a metallic powder by high pressure torsion can be properly modeled and the results can be used for material processing.


Introduction
High-pressure torsion (HPT) is an effective severe plastic deformation (SPD) process for producing materials with exceptionally small grain sizes ( ∼ 100 nm), and also for consolidating powder materials [1]. Grain refinement also occurs when powders are processed by HPT.
The majority of HPT studies are conducted on solid materials. However, there are obvious advantages of the HPT-assisted powder compaction for fabricating new materials. This technology not only provides efficient powder consolidation-even for hard-to-deform powders-but it also promotes the formation of strong contacts between particles, causes intensive mass transfer in materials that deform, activates phase transformations and chemical reactions [2]. These advantages recently generated a number of studies on HPT of powder materials (see, e.g. [3,4]).
The strain and porosity distributions in the volume of the specimen are very important for powder consolidation because these parameters determine the properties of the obtained materials. So far, this topic is quite scarcely studied. Experimental results on HPT consolidated water-atomized commercially pure iron powder were presented in Refs. [5,6]. In [5] it was reported that the porosity of the iron metallic powder saturates at about 3% under uniaxial compression. Torsion under axial pressure leads to further porosity reduction, up to tenths of a percent. In [6] it was shown that there is considerable strain heterogeneity in the axial direction of an HPT consolidated powder.
Mathematical modeling, employing, in particular, the finite-element technique, can provide useful information on HPT processing. Strain-induced phase transformations under compression and torsion in rotational diamond anvils were simulated using a finite-element approach in [7,8]. Significant heterogeneity was found and strong shear strain localization near the contact surface between the sample and the anvil was quantified. The effects of the elastic deformation of the anvil and the sample on the process were studied in [9]. The effect of friction, the thickness-to-diameter ratio of the specimen, and the rheology of the material on the stress-strain state of the sample during semi-constrained HPT were examined in [10,11]. The results showed that there is a tendency for flow localization during processing and this becomes more obvious during the processing of perfectly plastic or flow-softening materials or when processing samples having high thickness-to-diameter ratios. The analysis demonstrates the effects of the material condition, the disk aspect ratio, and the friction between the disks and the anvil walls. It has been shown that strain is non-uniform, not only in the radial but also in the axial direction of the sample. Extreme strain heterogeneity was reported in [12] where even a dead metal zone was found.
The elastic-plastic and rigid-plastic models are usually used for stress-strain state calculation in HPT. For example, the dislocation cell-based rigid-plastic model presented in [13] was used in [14]. A micro-macromechanical material model was employed in [15] for HPT simulations. None of the above-mentioned models allow investigating the evolution of the porosity during HPT of a powder material. The model in [7,8] takes into account the phase transformation with volume reduction, but the physical mechanism of volume change during the phase transformation is different from the one during shear of powders under pressure.
In HPT powder compaction, one has to take care about fully constraining the powder sample [1]. Unfortunately, all the above studies were done in unconstrained or semi-constrained conditions.
In this paper, we model the HPT of already precompacted powder specimens under fully constrained conditions; see previous publications of our group in [5,6]. For the simulation of the porosity distribution, we employ the model of [16]. Its main feature is that it takes into account the reduction of voids as well as the formation of new ones. The results show that by increasing the rotation angle of the anvil, the frontier of a high density area is spreading throughout the specimen. It starts from the edge of the rotating piston and moves in both radial and axial directions in the sample volume. We borrowed two model parameters from [16] for modeling the experimental results obtained in our group [5,6]. Good correspondence between theory and experiments was obtained concerning the strain and porosity distributions in the sample.

Modeling conditions
Several rheological models can be used in the finiteelement method for the simulation of plastic deformation of powder materials (see a review in [17]. The most commonly used are the modified Drucker-Prager model (see e.g. [18]) and the Shima-Oyane model [19] that are available in the ABAQUS and DEFORM software. In the present work, we apply the Beygelzimer model (see [16,20,21]), which combines the features of the Drucker-Prager model and the Shima-Oyane model. This model takes into account the reduction of voids as well as the formation of new ones. The basic equations of the model are the following (see e.g. [16]): where σ ik andė ik are the components of the stress and strain rate tensors, σ = 1/3σ ik δ ik andė =ė ik δ ik are the hydrostatic stress and the volume strain rate, respectively, while denote the square root of the second invariant of the deviatoric stress and strain rate tensors, respectively; θ is the porosity; k 0 is the generalized cohesion parameter and α is the generalized internal friction parameter from the Drucker-Prager model; a, m, n are the generalized parameters from the Shima-Oyane model. The technique to determine the parameters of the Beygelzimer model is presented in [16,21]. According the Beygelzimer model, the parameter α describes the structural ability of the material elements to accommodate each other. In the case when a complete accommodation of the elements is possible: α = 0. The value of α increases with the increase of difficulties to accommodate the plastic deformation. That is, the smaller the capacity of strain accommodation is, the higher the α is.
According to Zhao et al. [5], the porosity of precompacted pure iron powder was close to 3% before starting the rotation of the anvil. Torsion under pressure leads to a further reduction of the porosity. Therefore, the relative change in volume cannot be more than 0.03. Such a low porosity has a very small effect on the flow stress and the stress-strain state if there is no instability. Instabilities were not examined in our study. For this reason, the calculation of the stress-strain state of the sample was done using the approximation of a rigid-plastic incompressible body. The porosity was determined from the kinetic equation of the stress-strain state.
For materials with a low porosity and for compact materials, the value of the α is in the order of 0.01-0.001, a is in the order of 0.1-0.01, and k 0 = √ 2/3σ s , where σ s is the flow stress of the basic material (see e.g. [16]). Considering this, when θ 1, one can obtain from Equations (1) to (4), in the first approximation, the following equation system: whereė M = √ 2/3ġ is the strain rate and e M = t 0ė M dt is the true strain (von Mises strain).
The kinetic equation (6) is similar to the Shima-Oyane model [22] for the case when θ 1, however, it differs from it by the term α √ 3/2 on the right-hand side. This term is very important for HPT modeling because it take into account the formation of new pores during large shear strain.
According to Equations (5) and (7), one can find the stress-strain state of the specimen using the von Mises model for rigid-plastic incompressible materials (see, e.g. [23]). After that, the porosity of the material can be obtained by solving the kinetic equation (Equation (6)). This is, however, not possible when there is a loss of stability of the material because of the formation of shear bands (see e.g. [16,21]). We do not consider these regimes in this work.
The above-presented problem was solved by the finite element methods using the commercial software package DEFORM-2D/3D, V11.0 [24]. Figure 1 shows the schematic geometry of the constrained HPT process.
We have carried out a simulation of the HPT consolidation process according to the experimental conditions of our recent experiments published in [5,6]. The geometrical parameters of the simulations were the following: R = 10mm, H = 3 mm. According to Beigelzimer et al. [16], the α and a parameters for pure iron can be taken as α = 2.7 × 10 −3 and a = 2.7 × 10 −2 . The stress-strain curve obtained in our experiments [6] during the compaction of pure iron powder was built-in into the FE code ( Figure 2). The boundary conditions on the anvils were considered as perfect adherence of the powder to the flat contact surfaces, and the Trescatype friction law τ f = μσ s was used on the cylindrical surfaces (where μ is the friction coefficient; the value of μ was 0.25). Additionally, we suppose that there is no outward material flow during constrained HPT (see Figure 1). For this reason, the height of the specimen is not reduced during HPT. Under these conditions, a 2D axisymmetric torsion configuration could be used for the numerical calculations in the DEFORM software. The number of initial meshes of the workpiece was 10,000, which was sufficient to explore the local deformation behavior. There was no need to update the mesh in the DEFORM software because the mesh remained undistorted during the use of the 2D axisymmetric torsion model in DEFORM. The upper anvil was rotated at an angular speed of 0.1 rad/s up to two turns under a constant pressure of p = 1.1 GPa.
The porosity distribution was calculated using Equation (6) in 1000 nodes of the uniformly distributed square mesh. The results were examined in the axial crosssection of the specimen. The average porosity of the sample was calculated by the formula:

Results and discussion
The true strain distribution in the axial cross-section of the sample is shown in Figure 3. This figure shows that the HPT powder compaction process presents strong strain heterogeneities, not only in the radial direction, but also in the axial one. This result is in good agreement with our experiments [6]. For quantitative comparison, the calculated shear strain distribution at a distance of 9.61 mm from the axis of the sample after a rotation angle of π/2 is compared with the experimental one (see Figure 4). The discrepancy between experimental and theoretical strain near the top of the specimen (when z is small) may be due to (i) a large experimental error in the strain in this area (see Ref. [6]) and (ii) the slippage at the contact surface of the anvil. In order to examine the movement of the compacted powder zone, we take the true strain of 1.0 as a nominal value to define the border of a zone with large plastic deformation. Assuming proportional radial dependence of the shear strain, the true strain is defined for HPT by where β is the rotation angle of the anvil and r is the distance from the z-axis (see Figure 1). Figure 5 shows that this border moves in both the radial and axial directions with increasing rotation angle. Thus, the borderline depends on both r and z, meaning that Equation (9) is not valid for powder processing by HPT. It is not even valid for the r-dependence on the top surface of the disk where the e M = 1 borderline should be at the radial positions of 6.61, 3.31, and 0.83 mm for the rotation angles of π/4,  π/2, and 2π, while in the simulation it is at 5.5, 3.0, and 1.0 mm, respectively (see data in Figure 5, at z = 0). The model predicts a decrease in the average porosity with saturation at a level larger than zero. This is in quite good agreement with the experiment (see Figure 6). Note that the use of the non-modified kinetic equation of Ref. [22], which does not account for the appearance of porosity at large shear strain, leads to a very rapid decrease of the porosity. This is inconsistent with the experiments.
In the following, it will be shown that the movement of the border in the axial direction follows from the extremum principles of the theory of plasticity. Due to a small change in volumetric strain in comparison with deviatoric strain, it is possible to use the principle of minimum properties of an actual velocity field for incompressible material (see, e.g. [23]). According to this principle, the actual velocity field belongs to the minimum of the functional of the power dissipation W( V * ) of all admissible fields, which can be written for the present case as follows: Here the first integral is taken over the volume of the deformable body and the second is on the surface S f of its contact with the tool where there is friction between the deformable material and the tool. V * is a kinematically admissible velocity field, V * τ is its projection onto the surface S f , τ f is the friction stress andė * M is the von Mises strain rate corresponding to the kinematically admissible velocity field. The kinematically admissible velocity field must satisfy the boundary conditions for the velocity and also the incompressibility condition defined by div V * = 0 (see, e.g. [23]). Therefore, it can be taken in the form: where n ϕ is the unit vector in the direction of the tangential ϕ axis of the cylindrical coordinate system (r, ϕ, z) (see Figure 1), ω is the rotation rate of the upper anvil and (z) is an arbitrary smooth function of z, satisfying the boundary conditions: One can verify by direct calculation that the velocity field defined by Equation (11) satisfies the condition of incompressibility div V * = 0. The boundary conditions in Equation (12) express the condition that the material is sticking to the flat upper surface of the anvil and that there is only a tangential velocity component on the cylindrical surface at r = R. When (z) = z/H, we obtain the velocity field usually used for HPT. We generalize the expression of (z) by defining where k is a parameter. In the following, we will find out the value of k using the principle of minimum properties of an actual velocity field described above. For this purpose, we introduce the following dimensionless variables:z Here σ s0 is the flow stress of the material at the start of the rotation of the anvil. One can then show that the strain rate for the kinematically admissible velocity field has the form:ė * The projection of the velocity to the surface S f is given by One can rewrite then Equation (10) by using Equations (14)- (16) and the Tresca-type friction law At the start of the rotation of the anvil, the flow stress of the material is σ s0 , thereforeσ s = 1 throughout the volume of the sample (see Figure 7(a)).
One can get from Equation (17) for the start of the rotation of the anvil: By increasing k, the right-hand side of this ratio decreases monotonically, approaching to R/3 √ 3H when k → ∞. The dependence of the function (z) =z k on z when k 1 (fork = 100) is shown in Figure 8. It is clear from Figure 8 that only a thin material layer with a thickness of z 1 adjacent to the upper anvil is mainly deformed at the start of the rotation of the upper anvil. Plastic deformation leads to a nonuniform strain hardened layer (Figure 7(b)). We represent it as a homogeneous layer with an average flow stress of σ s1 (Figure 7(c)). The integrals over z in Equation (17) were done in two steps; over a segment z, then over a segment (1 − z). In this case, one can obtain from  Equation (17): For the case when z 1, this relation reads: The minimum of W is defined by the condition: When hardening is weak,σ s1 → 1, thus k → ∞. In this case, the velocity field becomes very heterogeneous and the deformation remains localized in the top thin layer. If the material hardens rapidly, the value of k is reducing according to Equation (22). For example, when σ s1 = 1.1 and z = 0.1 (the values of other parameters in Equation (22) are given above), one can get from Equation (22): k = 5.2. For this case, the (z) =z 5.2 function is displayed in Figure 8. One can see that for this value of k, the non-uniform torsion covers almost half of the height of the sample. Thus, the formation of a thin layer of the hardened material leads to the spread of the plastic deformation along the axial direction of the sample. Moreover, the strain rate in the hardened layer adjacent to the rotating anvil is significantly reduced.
As can be seen in Figure 2, the flow stress of the material is saturated at large strains. Therefore, at large strains, the flow stress becomes nearly uniform within the sample volume. This effect begins in the upper layer, where the strain is the largest. It is well known that low hardening rate can lead to strain localization. For these reasons, both at small and at high rotation angles, the deformation is localized (see Figure 3). One can see that only a thin material layer adjacent to the upper anvil is strongly deformed at the start of the rotation, followed by plastic deformation in the whole volume, finally, at large rotation angles, the deformation localizes again in the upper layer of the specimen.
The above-discussed development of the plastic deformation in constrained HPT testing follows from the principle of minimum power dissipation. Since this analysis is based on the von Mises model, the obtained results are valid not only for low porous materials but for bulk materials too.

Conclusions
By employing the Beygelzimer model of a structurally inhomogeneous porous body (see [16,20,21]) in the commercial DEFORM-2D/3D V11.0 finite element software, it is shown in this paper that the strain under constrained HPT is non-uniform in both the radial and axial directions. For initially homogeneous samples, the deformation is concentrated in a thin layer adjacent to the rotating anvil when rotation starts. If the material does not display strain hardening, the deformation remains localized in a thin layer during further rotation. In contrast, if the material is strain hardening, the plastic deformation penetrates into the sample progressively as a function of the rotation of the anvil. The border of the plastic zone with large strain moves in both the radial and axial directions by increasing the angle of the anvil rotation. The obtained effects are valid for both powder and bulk materials. These effects are associated with the friction on the lateral surface of the sample with the anvil and follow from the principle of minimum power dissipation.
The mathematical simulation predicts a decrease in the average porosity with saturation at a level more than zero. This is in quite good agreement with our experiments recently published in [6].