Calculation of the Plastic Metal Flow in the Cold Extrusion Technology by Finite Volume Method

This paper presents an application of the finite volume method to ideal plastic metal flow in extrusion technology. Governing equations for the mass and momentum balance are used in an integral form. Solution domains in the cases considered are discretized with a Cartesian numerical mesh with computational points placed at the center of each control volume. After discretization of the governing equations, the resulting system of nonlinear algebraic equations is solved by an iterative procedure, using a segregated algorithm approach. Resulting stress fields are obtained from the Levy-Mises equations. The experimental results and numerical calculations are in good agreement.


INTRODUCTION
Extrusion processes in the metal forming by plastic deformation belong to the widespread production technology. One of the important parameters of this process is the working pressure, i.e. deformation force that determines the shape and dimensions of the extrusion tool and characteristics of the machine on which the technology will be performed. Also, a more complete knowledge of metal flow in these processes allows better design of tooling, especially in the extrusion of rods through the dies with complex cross sections.
In order to obtain a high-production rate with acceptable quality, many process parameters must be controlled. The plastic flow is conditioned by the properties of the extruded material, the geometry of the tool and the condition of the workpiece-tool contact surfaces (friction). The paper discusses the flow of ideally plastic material at plain-strain extrusion. Forward extrusion is a forming process in which a workpieces pushed through a die whose exit diameter is smaller than that of the work piece.
Thirty years ago, Demirdžić et al. [6,7] proposed the first application of the cell-centered finite volume method, in its modern form, to solid mechanics. This original method was applied to the simulation of thermal deformations in welded workpieces on a 2-D structured quadrilateral grid. Small strains and rotations were assumed and the material behavior was described by the so called Duhamel-Neumann form of Hooke's law [6].
This cell-centered approach takes its name from the dependent variable, in this case velocities, residing at the cell centers (control volume centroids). Calculations performed in this paper were done with Comet software (Version 1.052) [8].

MATHEMATICAL MODEL 2.1 Governing Equations
The theory of plastic flow in the problems of bulk metal forming uses point velocities as dependent variables instead of displacements. The approach to problem solving according to this theory is analogous to the approach in fluid mechanics. The cold extrusion processes, which are considered in this paper, assuming a quasi-stationary (steadystate) flow, are also solved according to this theory and must satisfy the following basic conservation equations: Mass balance equation (continuity equation) [7]: Momentum balance equation [7]: Eq. (1) and (2) are valid for an arbitrary part of a continuum of volume V bounded by the surface S. Here t is the time, ρ is the density, u j is the velocity vector, σ ij is the stress tensor and f bi is the resultant body force.

Constitutive Relations
According to the theory of plastic flow, the deviators of stress tensors and strain rate tensors are coaxial. The components of stress tensor are determined by the Levy-Mises equations [1]: where σ is the mean stress (pressure), while the viscosity µ can be calculated as [1]: Whence it follows that the Levy-Mises equations hold for 0 ≠ ε otherwise the material behaves like a rigid body.
The effective stress σ and effective strain-rate ε are defined by equations [1]: where: ij σ ′ is deviatory stress and ij ε strain-rate tensor with components [1]: The viscosity µ is in the general case a function of the degree and rate of deformation and temperature, i.e.
Using the assumption of ideal plasticity and the fact that cold forming is considered, viscosity becomes only a function of the strain rate gradient. The initial viscosity µ p at which the plastic flow begins is determined by the uniaxial tensile experiment (σ z = σ y = 0, ).
In this case, Eq. (5) and Eq. (6) become: , so the viscosity at the beginning of plastic flow in the direction of x-axis will be [1]: where σ Y is the yield stress, while the strain rate , h u x = ε is equal to the ratio of the ram speed u, m/s, and workpiece high h. The value of the viscosity calculated according to Eq. (8) is also the upper limit that the viscosity can reach in the calculation. Areas in the deformed material where the viscosity is above this limiting value were considered in the calculation as rigid (stagnation) zones.
By entering the constitutive relation (3) into the momentum equation (2), a closed system of equations is obtained, which in the general 3D case consists of four equations with four unknown functions (u, v, w and σ) of spatial coordinates and, for nonstationary cases, time.

The Finite Volume Discretization
In order to solve the defined system of equations, it is necessary to discretize time, spatial domain and the equations themselves. The general case of discretization of Eq. (1) and Eq. (2) is described in [7]. Two -dimensional discretization is considered in this paper.
With the assumption of incompressibility of the extruded material, as well as the previously introduced assumption of stationarity of the process, the first terms on the left side of Eq. (1) and Eq. (2) become zero. Also, the effect of body forces, that are present only in so-called high-speed metal forming, is neglected.
The two-dimensional spatial domain (the domain of solution) is sub-divided into a finite number (N) of contiguous control volumes (CV) or cells of volume V bounded by the surface S consisting of a number of cell faces of the area Sk. The solution domain is discretized by a rectangular numerical grid, example on Fig. 1.   Fig. 3, [9]. The calculation points are located in the centers of gravity of the control volumes, while boundary nodes, needed for the specification of boundary conditions, reside at the center of boundary cell faces.
To illustrate the method of discretization of Eq. (2), the equation of balance of forces acting on one control volume, Fig. 4, will be set only for x direction.
Based on Eq. (2) for the stationary case and in the absence of volumetric forces, the equilibrium equation for one CV can be written as [9]: Where k = e, n, w, s is the counter of the sides of the control volume. The last term in Eq. (9) appears only for axisymmetric cases of the stress state, in the equilibrium equation for the x direction. Based on Eq. (9), the forces on individual sides S k of the observed control volume are given by the following expressions [ The terms on the left side of Eq. (11) represent the convective change, while the terms on the right side are the acting surface forces on the control volume.
Convective terms are discretized by applying a mixture of formulas for upstream differentiation of the first order of accuracy (implicit part) and the difference of symmetric central differentiation of the second order of accuracy and upstream differentiation (explicit part), [7].
The values of the variables between two adjacent calculation points (on the sides of the CV) are obtained by linear interpolation.
The terms on the right hand side in Eq. (11), so-called diffusion terms are discretized as follows [9]: The coefficients for creating a system of algebraic equations in order to solve unknown velocity components in the x direction have the form Each of four coefficients represent the sum of the implicit part of the diffusion terms and the implicit part of the convective terms for the observed CV, where k m  (k = e, n, w, s) is the mass flus through certain CV surface.
The so-called source terms (b) of the system of algebraic equations are the sum of the explicit parts of the convective and diffusion terms and the terms of the forces from the mean normal stress σ (pressure).
The final appearance of the discretized equation for one CV for calculation of the velocity component in the x direction is [9]: ( , , , ) The discretized equation for the y direction is analogous to Eq. (14), where the coefficients a K (K = E, N, W, S) and a P remain the same as for the x direction. The difference is only in the source term b.

Calculation of Pressure
It is easy to notice that in the discretized equation of momentum there is an unknown mean normal stress (pressure), and that the equation of momentum (Eq. (2)) is not independent of the continuity equation. However, the pressure does not figure directly in the continuity equation, but it is used as an additional condition that must be satisfied when calculating the velocity field. Therefore, several procedures have been developed to connect the equations of momentum and continuity by which the pressure is determined. For calculation in this paper, the SIMPLE algorithm [7] was used, which results in a linear algebraic equation for pressure correction that also has the form of Eq. (14).

Boundary Condition
To start the calculation, all dependent variables (velocity components) have to be set to their initial values. The expressions for the evaluation of forces are valid for all interior cell faces which are shared by the two adjacent cells. Boundary cell faces are related to only one cell and require special treatment that depends on the type of the boundary conditions. Boundary conditions can be classified as being of Dirichlet (e.g. 'inlets', where all dependent variables are usually known), or Neumann type where values of dependent variable gradients are available. If the problem considered is symmetric in such a way that one boundary can be taken for a symmetry plane, then that is possible by looking at a near-boundary cell and its mirror image on the other side of the symmetry plane.

Solution Procedure
Writing the equations of form (14) for each velocity components and for each CV, a system of 3×N (N -number of CV-a) of nonlinear, coupled equations is obtained, which is solved by the following iterative procedure [7].
The equations are temporarily linearized and decoupled assuming that the coefficients and the source terms are known (calculated based on the values from the previous iteration) thus obtaining 3 systems of N linear algebraic equations. The coefficient matrices for each subsystem of linear equations will be sparse, symmetric, and diagonally dominant. These subsystems of the equations are solved one after the other by the conjugate gradient method. After that, the coefficients and the source terms are renewed and the procedure is repeated until the desired accuracy is achieved.

NUMERICAL SIMULATION OF EXTRUSION PROCESS
To verify the presented numerical procedure, two cases of the plain-strain extrusion (forward and backward) will be considered. Experiments for these cases were done by Segal [10]. In numerical simulation these processes are considered as steady-state.
The first case is forward extrusion under ideal contact conditions (contact tangential stress τ k = 0) with a cross section reduction ratio of A 0 /A 1 = 2. The initial width of the piece is 40 mm, while the width at the output cross section is 20 mm. The discretization of the spatial domain was done with a rectangular numerical grid with 430 CVs, Fig. 5. The effective strain-rate distribution over solution domain is given in Fig. 7. The areas which are not covered by contours represent the stagnation (dead) zones. Fig. 8 shows the fields of the velocity components in the x-direction obtained experimentally by the Moire method Symmetry plane [10], and by the finite volume method. There is a very good accordance between the calculated and experimental results.  Another example relates to the backward extrusion of the same material and process parameters, under a plain-strain and no friction conditions, Fig. 9. The initial width of the prism is b 0 = 80 mm, while extrusion is done by using a b 1 = 44 mm wide punch. The discretization of the spatial domain was performed with 450 CVs. Fig. 10 and Fig. 11 show the fields of the velocity components in the x and y-direction obtained experimentally, by the Moire method [10], and by the FVM. Accordance between the results is very good.  The assumption of the existence of rigid (dead) zones during the deformation process is one of the common assumptions in solving the problem of the theory of plasticity by using the concept of a rigid-plastic body (for example slipline field method, upper bound method, thin section method), [1,2,11,13]. Rigid zones are most often assumed in the corners of extrusion dies (see Fig. 6, Fig. 8) or under punch (see Fig. 10, Fig. 11). The shape of these zones obtained by measuring the velocity components experimentally as well as by the FVM, indicates the accuracy of these assumptions. The height of these zones increases with increasing the contact friction between the workpiece and the die.
In the flow approach, the material is assumed to behave like a non-Newtonian viscous or visco-plastic fluid and the numerical solution is obtained by using an Eulerian reference frame (stationary mesh).

CONCLUSION
The finite volume method can be successfully applied in the analysis of plastic flow of metals in metal forming technologies. In this paper, the method was used for the analysis of the flow of ideally plastic material and the orthogonal geometry of tooling. Unlike approximate methods of plasticity theory, FVM gives the complete picture of the considered flow, i.e. arrangement of velocity components throughout the spatial domain, from which other relevant variables can be easily calculated: strain rates tensor components, stress tensor components, effective stress distribution, etc. Compared to other numerical methods, the advantage of FVM is in more efficient use of computer resources due to the structure of the matrices of a system of algebraic equations and solving systems of linear equations separately.
The further application of FVM in the analysis of plastic flows of metals could be: inclusion of arbitrary geometric shapes of solution domain (including 3D domains), calculation the effects of material hardening, inclusion of friction on contact boundaries and analysis of hot extrusion processes.
FVM can be applied to the above cases by generalizing the presented discretization procedure and including the energy balance equation [8]. In that case, a new subsystem of linear algebraic equations is obtained with unknown temperatures in the CVs centers. Some results of FVM application for non-orthogonal solution domains and nonzero friction boundary conditions are presented in [8].

Notice
The paper was presented at MOTSP 2020 -International Conference Management of Technology -Step to Sustainable Production, which took place from 30 th September -2 nd October 2020 in Bol, island Brač (Croatia). The paper is not and will not be published anywhere else.