Finite Volume Analysis with the MacCormack Method Applied to Metal Flow in Forward Extrusion

Computational numerical simulation has been largely applied in the design and analysis of metal forming processes. Extrusion is one of the main forming processes largely applied in the manufacturing of metallic products or parts. Historically, the Finite Element Method (FEM) has been applied for decades in metal extrusion analysis. However, recently in the academy, there is a trend to use the Finite Volume Method (FVM), because literature suggests that metal flow by extrusion can be analyzed by the flow formulation. Thus, metal flow can be modelled as an incompressible viscous fluid. The MacCormack Method is commonly used to simulate compressible fluid flow by the FVM. However, metal extrusion does not present state equations to calculate the pressure, and therefore, a velocity-pressure coupling method is necessary to obtain consistent velocity and pressure fields. This work proposes a new numerical scheme to obtain information about metal flow in the extrusion process along the steady state. The governing equations were discretized by FVM, using the Explicit MacCormack Method to structured and collocated mesh. The SIMPLE Method was applied to attain pressure-velocity coupling. This new numerical scheme was applied to analyze forward extrusion of lead. The metal extrusion velocity fields were calculated with a fast convergence and presented a good agreement with analytical and experimental results obtained from literature.


Introduction
In the academic and industrial environments, numerical simulation is an unquestionable reality. This can be seen by a great number of scientific papers published in journals and conferences that discuss the theme [4,8,10,15,17]. In addition, many industries use commercial codes or codes made by their own engineers to analyze projects [5,8,15,17]. In the metal forming area, it is not different. There are many scientific papers published to evaluate metal forming phenomena by numerical simulation [1,6,7,9,11,12,14,18]. This is because computer has steadily increased the processing power since the 1950's [4,8,15,17]. Simultaneously, the numerical methods applied to solve mathematic equations governing physics phenomena of metal forming, had to be adapted to specific situations [4,8,15,17]. On this way, one can highlight the importance of Finite Element Method (FEM) to simulate metal forming process [8,15,17].
The FEM was created in the early 1950's linked with elasticity structural problems. In the next decade, the FEM was applied in metal forming. But at that time, it was utilized solid formulation, considering elastoplastic materials and elasto-viscous-plastic materials. However, this formulation was inadequate to analyze the nonlinearities (geometric and material) which appeared from these problems [8,15,17].
In 1970s, it was developed the flow formulation. In metal forming processes, the flow formulation considers the metal similar to an incompressible viscous fluid [15,17], which was and still is analyzed numerically by Finite Volume Method (FVM) [4]. The flow formulation can be assumed because metal forming processes are isochoric, i.e., with constant volume [15,17]. The flow formulation considers that materials are rigid-elastic-plastic and rigid-viscous-plastic. Thus, the elastic deformations are neglected because the high value of plastic deformations in metal forming processes [8,15,17].
Nowadays, the FEM is still a preferential choice to analyze metal forming problems [8,15,17]. However, from the 1990sthere was an increase of interest in using FVM in metal forming problems [6,7,9]. This happened because FVM is a method which is well tested in fluid dynamics and this knowledge would help to understand metal forming phenomena and create numerical codes which can involve analysis of both solid and fluid or interactions between them [4]. The FVM was created in the 1970s for solving the limitations of Finite Difference Method (FDM) [4,10] and since then has been the numerical method most used in fluid dynamics [4,10]. Increase in utilization of FVM came from its main feature, that is, to satisfy the conservations laws in global and discrete levels [4,10]. Thus, it is linked with more consistent physical interpretation. Knowing that, FEM and FDM do not satisfy conservations laws in discrete level because they do not work with control volumes but with mesh points [4]. Then, from the 1990s, a lot of works appeared in literature to numerically simulate metal forming process by FVM and flow formulation [6,7,9,11,12,14,18].
In the numerical point of view, fluid flow can be classified as compressible (when the density is a strong function of pressure) or incompressible (when the density is a constant or is an only function of temperature) [4]. Therefore, to analyze an incompressible fluid flow, the parameters used are velocity, temperature and pressure fields [4,10,13]. In metal forming flow, without significant variation of density, it is necessary to determine these same fields [15,17].
Considering incompressible fluid flow and metal flow, both cases are governed by conservation law equations (mass, momentum and energy), those are second order partial differential equations (PDEs) and can be represented by Generic Transport Equation [2,4,10,13], which has four parts, being: transient term, convective term, diffusive term and source term [4,6]. Each one of these terms represents the physical behavior of the phenomena and these behaviors are a mix of elliptic, parabolic and hyperbolic behavior [4,10].
The Generic Transport Equation represents a system of partial differential equations. In the numerical solution, this system is transformed into a set of non-linear algebraic equations if an implicit approach is applied or into a set of non-linear equations if explicit approach is applied [4]. To solve these systems or set of equations, two ways can be applied: coupled or segregated solutions [4]. If the segregated solution is chosen, the system of each variable is solved in the sequential manner, but there is a problem with coupling among some variables, mainly velocity and pressure, because there is no equation to calculate the pressure field, by the fact is that the density has to be calculated by state equations [4,10].
Besides that, if collocated mesh is to be used, what is more frequently, it is raised the problem of "checkerboard" oscillations, because the computational point used to calculate the velocity and the pressure is the same [4,10]. To solve these problems, it is necessary to apply a velocity-pressure coupling method, such as, SIMPLE, PISO, etc., and a method of interpolation to move the velocity components to different positions from the pressure ones [4,10].
In order to close the set of equations needed to solve the metal flow problem, the constitutive equation must be added, which, of course, is different in the metal flow and in the incompressible fluid flow. The governing equations and the constitutive equations together are necessary and sufficient to consistently apply a numerical method [4].
There are many approaches to generate approximation equations by FVM [4,10], and all methods should satisfy Lax's Equivalence Theorem to assure convergence for the numerical method applied to solve PDEs [10,11,14].
The MacCormack method is traditionally used in compressible fluid flow in aerodynamics analysis, because it smooths discontinuity peaks produced by shock waves [10]. The MacCormack method is a two-step method of the Lax-Wendroff scheme and its main features are to have second order accuracy in time and space [10].
The MacCormack Method also presents good results in the solution of non-linearities of Partial Differential Equations [10]. As metal forming flow is governed by Conservations Differential Equations, this work proposes to use the MacCormack method to solve the governing equations of metal flow in forward extrusion, considering the velocities discontinuities at the entrance and at the exit of the deformation zone [15].

Governing Equations
The mathematical modelling in metal flow of forward extrusion can be treated by FVM using flow formulation [1,6,11,14,15,18]. According [15] metal flow is considered similar to fluid viscous flow. Reference [10] suggests fluid viscous flow can be mathematically represented by mass, momentum and energy conservation equations (1), (2) and (3) respectively.
Whereρ is density, v is velocity vector, b is body forces, σ is Cauchy stress tensor, q  is heat flux vector, e is specific energy, ρ is heat generate and t is time [2,10]. The equations (1) to (3), can be written in the same form, called the generic transport equation. This equation governs the transport of the φ property inside of matter and it is given by [6]: Where φ is the property transported and φ Γ is a diffusive term [6]. The first right side term is the transient term and represents the accumulation of property inside the control of volume and commands the time variation. The second right side term is the convective term and represents the property flux leaving the control volume due velocity and it has a Universal Journal of Mechanical Engineering 5(1): 1-8, 2017 3 hyperbolic/parabolic nature. The first left side term is the diffusive term and represents the property flux leaving the control volume due to molecular diffusion and it has elliptic nature. The last left side term represents the source term, i.e. inside the source term are all terms to represent whether the property is created and destroyed or other terms are not convective, diffusive and transient [4,1013].
Considering (1), (2) and (3) in cylindrical coordinates for an axisymmetric case, these equations can be written in matrix structure in the following form [10,13,14]: being ρ is radial coordinate, z is axial coordinate, Q, F r , F z e S are flux vectors that assume the following format [10,14]: With σ ρρ , σ θθ and σ zz normal components and σ ρz the shear component of Cauchy stress tensor, T is the temperature, c is the specific heat, ρ q  and z q  are flux heat vectors, v ρ and v z are velocity vector components, σ the effective stress and ε  is the effective strain rate [10,15,17].

Plasticity Constitutive Equations
Plasticity constitutive equations are mathematic expressions that describe the material behaviour under mechanical and/or thermal loading by many variables: static (given by stresses), kinematic (given by displacements, deformations and velocities) and thermal (given by flux thermal and temperature) [13]. Metal constitutive equation for axisymmetric plastic flow assumes the material to be incompressible, isotropic and rigid-viscous-plastic [14,15].
The relation between the stress tensor and the strain rate tensor can be written in the form of (10) [6]: (10) where σ is the Cauchy stress tensor, ε  strain rate tensor, σ m hydrostatic pressure and η metal viscosity. Hydrostatic pressure inside the solid is calculated by (11) [13]: where, σ ρρ , σ θθ and σ zz are principal stresses of the stress tensor. Metal viscosity η is defined by the Associated Plastic Flow Rule given by (12) [15,17]: being σ ′ the deviatoric tensor and λ  the plastic multiplier. Metal viscosity η can be calculated by (13) [15,17]: where σ is the effective stress and ε  is the effective strain rate. The effective strain rate must not assume values less than 10 -3 to avoid numerical instability [6].

Finite Volume Method
Considering (5) and integrating it over the control volume given by Figure 1A, using Gauss' Theorem, results in (14) [10,14]: Where S is the outward surface vector and V mn is the control volume area.

MacCormack Method
Explicit MacCormack scheme on variable time was applied in (14). Considering the outward surface vector on the control volume given by Figure 1b, (15) represents the predictor step; (16) represents the corrector step and (17) represents the current step [11,10,14]. In this study, the MacCormack scheme is a pseudo-transient process, where ∆t is a virtual time increment to obtain the final converged solution [10].  (17) such, t is current time, (t+1) is next time step, ∆t is a virtual time step, is corrector step and 1 + t mn Q is current step. The flux vector F, inside (15) and (16), must be discretized in correct way to ensure the main feature of MacCormack Method: second order accuracy in time and space [10].

Velocity-pressure Coupling
The SIMPLE Method was applied to obtain velocity-pressure coupling in the numerical scheme developed in this work [4,10]. The SIMPLE Method makes the coupling by corrections in both variables, i.e., velocity and pressure. To attain this, the method utilizes (18) Where the variables p and v are current values of pressure and velocity vector, respectively, and the variables p 0 and v 0 are the estimated pressure and velocity vectors and the variables p´ and v´ are corrections of each variable. According to [4], the corrections equations in discretized volumes are given by: Where lowercase subscripts represent the face position of control volumes and upper case subscripts represent the centre position of control volumes, according to Figure 2 [16]. The term A is defined as the quotient between fictitious time step and density, according to [10]. The velocities v ρ and v z must be considered on the faces of control volumes, i.e., the velocities must be dislocated from control volume centre to the control volume faces in order to have the ideal coupling between velocity and pressure fields. To this process, it was used in the Momentum Interpolation Method Universal Journal of Mechanical Engineering 5(1): 1-8, 2017 5 [3]. Reference [4] suggests the correction equation (24) to pressure or the Poisson equation that calculates p´: where: being ρ n the radius on n face and ρ s is radius on s face of control volume.

Boundary Conditions
The boundary conditions were imposed over the control volume faces, using the scheme called "Ghost Volume". In these ghost volumes boundary conditions of "Diρichlet Kind" and "Neumann kind" can be applied. According to [4,6] the boundary conditions applied for metal flow in forward extrusion, analysis should be imposed in this way, being control volumes walls represented by Figure 3: where t nt is friction stress, μ is Coulomb coefficient and σ n is contact normal stress.

Case Description
The metal flow analyzed in this work refers to the forward extrusion of lead with die geometry described in reference [6] and shown in Figure 3, where inlet radius, outlet radius and die semi-angle of 15 o , define a cross-section reduction equal to 25%. Figure 3 also shows the mesh used in the simulations, with 1463 quadrilateral volumes. The extruded material was considered to be rigid-perfect-plastic, and the Coulomb friction coefficient was assumed in a range from zero to 0.5 and applied in the conical die region. Other parameters of the forward extrusion and material are presented in Table 1.

Discussion of Results
The numerical results obtained in this work from the  [5,8,15]. The numerical results from literature are presented in the paper [6]. Figure 4A shows a comparison between the results of radial analytical velocity (v ρ_analytic , generated by Upper Bound Method) and radial numerical velocity (v ρ_numeρic , generated by present numerical scheme). As can be seen, the values and contours are in good agreement. Figure 4B shows a comparison between the results of axial analytical velocity (v z_analytic , generated by Upper Bound Method) and axial numerical velocity (v z_analytic , generated by present numerical scheme). Again, the results of values and contours are in good agreement. Figure 5A shows a comparison between axial numerical velocity contours (V z ) of present work and axial numerical velocity contours (V z ) presented in [6]. It can be noticed that values and contours are in good agreement. In Figure 5B is showed a comparison between radial numerical velocity contours (V ρ ) of present work and radial numerical velocity contours (V ρ ) presented in [6]. Again, the radial numerical velocity results are in good agreement, but the contours have some differences which can be assigned by different ways that were applied in the Coulomb friction model. Because in the present work, Coulomb´s friction coefficient assumed the normal pressure at boundary, while in the paper by [6] it was assumed the hydrostatic pressure. Figure 6A presents pressure contours, where, it can be observed that there is a compressive region at the conical die region, close to the material-die interface. Figure 6B shows effective strain rate contours, where it can be seen that the biggest effective strain rate arise at conical region close to the material-die interface. Figure 6C shows the effective stress contours, where it can be observed rigid regions of material during forward extrusion process.  of present work and axial numerical velocity contour presented by [6]. (B) Comparison between radial numerical velocity contour (Vρ) of present work and radial numerical velocity contour presented by [6].

Conclusions
The present numerical scheme applied in metal flow of lead forward extrusion allowed the following conclusions:  Numerical results were calculated by the present finite volume numerical scheme without artificial viscosity;  The numerical scheme presented is explicit, then it is conditionally stable because it was used a time step equal 10 -15 to satisfy CFL (Courant-Friedrichs-Lewy) condition;  the axial velocity increased at the entrance of the conical die region, and the biggest axial velocity was achieved at the outlet die region;  at the conical region of the material-die interface occurred the biggest pressure and the biggest effective strain rate;  the numerical scheme convergence was achieved, for axial and radial velocity and pressure, with about 20,000 iterations;  flow formulation with both FVM and MacCormack Method was able to produce good results for metal flow of lead forward extrusion;  MacCormack Method besides to be applied in compressible flow analysis, as proved by the literature, can also be applied to model metal flow in forward extrusion.