Efficient and accurate approach for powder compaction problems

Efficient and accurate approach for powder compaction problems A. Pe ´rez-Foguet, A. Rodrı ´guez-Ferran, A. Huerta Abstract In this paper, a new approach for powder cold compaction simulations is presented. A density-dependent plastic model within the framework of finite strain multiplicative hyperelastoplasticity is used to describe the highly nonlinear material behaviour; the Coulomb dry friction model is used to capture friction effects at die-powder contact; and an Arbitrary Lagrangian–Eulerian (ALE) formulation is used to avoid the (usual) excessive distortion of Lagrangian meshes caused by large mass fluxes. Several representative examples, involving structured and unstructured meshes are simulated. The results obtained agree with the experimental data and other numerical results reported in the literature. It is shown that, contrary to other Lagrangian and adaptive h-remeshing approaches recently reported for this type of problems, the present approach verifies the mass conservation principle with very low relative errors (less than 1% in all ALE examples and exactly in the pure Lagrangian examples). Moreover, thanks to the use of an ALE formulation and in contrast with other simulations, the presented density distributions do not present spurious oscillations. Keywords Powder compaction, finite strain multiplicative plasticity, consistent tangent moduli, Arbitrary Lagrangian–Eulerian (ALE) formulation, density-dependent models, numerical differentiation 1 Introduction Cold compaction processes are a key ingredient in powder forming processes. They consist in the vertical compaction through the movement of a set of punches of a fine powder material at room temperature. The process transforms the loose powder into a compacted sample with a volume reduction (and therefore a density increase) of about 2–2.5 times. The design of these processes includes the definition of the initial dimensions of the sample and the movements of the punches that lead to compacted samples with uniform density distributions. In this context, efficient and reliable numerical simulations can play an important role as a complement of experimental tests. Two ingredients are crucial for the numerical modelling of powder compaction processes: the constitutive relationship and the kinematic formulation of the problem. Several constitutive models have been proposed, including microscopic models, flow formulations and solid mechanics models, such as elastic, plastic or viscoplastic models; see Oliver et al. [1] and Lewis and Khoei [2] for a general overview and references for each type of model. One of the most common approaches is the use of elastoplastic models based on porous or frictional materials. Here, plastic models expressed in terms of the relative density and the Kirchhoff stresses are considered. The constitutive relationship is formulated within the framework of isotropic finite strain multiplicative hyperelastoplasticity [3, 4]. This type of models has already been applied to powder compaction problems with some simplifications derived from the assumption of small elastic strains [1]. In this work, large elastic strains are included in the formulation [5]. It is shown that this does not represent any drawback from a modelling point of view. On the contrary, it allows to apply numerical techniques and material models developed for the general kinematic framework in a straightforward manner. Up to date, a common feature of powder compaction simulations with solid mechanics constitutive models is the use of a Lagrangian kinematic formulation. This approach has shown to be adequate for problems that do not exhibit large mass fluxes among different parts of the sample (i.e. homogeneous tests). But in practical problems, as those which appear in realistic design processes, the Lagrangian approach leads to too distorted meshes [2, 6] or violations of the boundary limits [1]. In order to solve these problems, different h-adaptive procedures have been presented recently [6, 7]. However, h-refinement is computationally expensive and information must be interpolated from the old mesh to the new mesh. For these reasons an Arbitrary Lagrangian–Eulerian (ALE) approach is preferred in this work. ALE formulations were first proposed for fluid problems with moving boundaries [8, 9]. Nowadays, ALE formulations for fluid problems are widely used in forming processes. On the other hand, the ALE formulation has been successfully employed in nonlinear solid mechanics [10–14]. The ALE formulation for Computational Mechanics 30 (2003) 220–234 Springer-Verlag 2003 DOI 10.1007/s00466-002-0381-4 220 Received: 20 March 2002 / Accepted: 15 October 2002 A. Pe ´rez-Foguet (&), A. Rodrı ´guez-Ferran, A. Huerta Depart. de Matema `tica Aplicada III, E.T.S. de Ingenieros de Caminos, Canales y Puertos, Universitat Polite `cnica de Catalunya, Jordi Girona 1, E-08034 Barcelona, Spain e-mail: agusti.[email protected] The partial financial support of the Ministerio de Ciencia y Tecnologı ´a (grant number DPI 2001-2204) is gratefully acknowledged. multiplicative hyperelastoplasticity recently presented by Rodrı ´guez-Farran et al. [14] is used here. Due to the high nonlinearity of powder compaction simulations, a common approach is the explicit treatment of both the integration of the constitutive equations and the solution of the equilibrium equation [2, 6, 15], although, implicit schemes have also been successfully applied [1, 7]. In fact, the development of robust, accurate and efficient schemes, both implicit and explicit, is still an active research topic, see for instance [16–19] among many others. Here a mixed approach is followed. The integration of the ALE constitutive equations is done by means of a fractional-step method [3]: the treatment of the Lagrangian phase is implicit and the scheme for the convective phase is explicit. The equilibrium equation is solved with an incremental-iterative approach, where only the Lagrangian phase is performed within the iterations of each load increment (the perturbation due to the integration of the convective phase is taken into account in the subsequent load increment). A key point for the efficient implicit solution of the equilibrium equation is the consistent linearization of the algorithm, which is given by the proper consistent tangent moduli [20]. The recently presented expression for density-dependent models by Pe ´rez-Foguet et al. [5] is used here. Specific tools for two of the main issues of implicit approaches for complex elastoplastic models are also considered within this work: the convergence of the plastic corrector at Gauss-point level is guaranteed by the use of globally convergent Newton–like schemes [21, 22], and numerical differentiation schemes are used for the computation of complex plastic equations derivatives following references [18, 23]. An outline of the paper follows. The main features of the constitutive equations and the numerical time–integration scheme are presented in Sect. 2. The consistent tangent moduli is presented in Sect. 3. After that, in Sect. 4, the proposed approach is applied to several representative powder compaction problems. The present results are compared with experimental data and other results of numerical simulations presented in the literature. Section 5 contains some concluding remarks. 2 Problem statement In this section, the basis of the proposed approach is briefly presented. First, the kinematics of the ALE formulation are reviewed. After that, the finite strain elliptic elastoplastic model is described, focusing in the yield function dependence on the density. Finally, the numerical time-integration algorithm is presented. 2.1 Kinematics Let RXRndim (ndim ¼2;3) be the material configuration of a continuum body with particles labelled by their initial position vector X2RX. In a Lagrangian setting, Xare used as the independent variables in the description of motion. The motion of the body is described by the one-parameter family of mappings ut:RX7! Rndim with t2½0;T. Rx¼utðRXÞis the spatial configuration of the body at time t, and x¼utðXÞ¼uðX;tÞ2Rxis the current position of the material particle X. The key ingredient of the ALE formulation is the referential configuration Rv, with grid (or reference) points v used as independent variables to describe body motion. This referential configuration Rvis mapped into the material and spatial configurations by Wand Urespectively, see figure 1. The initial position of a material particle is expressed on the reference domain as X¼Wðv;tÞ, and the current position as x¼Uðv;tÞ. The three mappings u,U and Ware related by u¼UW1. In an ALE setting, different displacement fields can be defined. Two of them have special interest, the particle displacement uand the mesh (spatial) displacement uU, uðX;tÞ¼xðX;tÞXand uUðv;tÞ¼xðv;tÞv: ð1Þ The particle velocity and the mesh velocity are respectively v¼ox otjXand vmesh ¼ox otjv;ð2Þ where, following standard notation, jmeans ‘‘holding  fixed’’. The link between material and mesh motion is provided by the convective velocity, c¼vvmesh. 2.2 Constitutive model Isotropic finite strain multiplicative plasticity is assumed. The particle deformation gradient, FðX;tÞ¼ou oXðX;tÞ;ð3Þ is locally decomposed into elastic and plastic parts as F¼FeFp. The thermodynamic state is defined by means of the elastic left Cauchy–Green tensor be¼FeFeT, where the superscript T means transpose. As usually in densitydependent plasticity, no plastic internal variables are considered [1, 5]. The Kirchhoff stress tensor, s, is given by the hyperelastic relationship Fig. 1. Domains, mappings and deformation gradients in the ALE description 221 s¼2dWe dbebe;ð4Þ where Weis the free-energy function per unit of undeformed volume. The Cauchy stress tensor is given by r¼s detðFÞ:ð5Þ The plastic response of the material is assumed isotropic and density-dependent. The dependence on the density is incorporated in the yield function fðs;gÞand the flow direction mðs;gÞthrough the relative density gðX;tÞof the particle Xat time t. The relative density is equal to the real density of the material divided by a reference value, which usually is the solid density of the compacted material. Associative plasticity is characterized by mðs;gÞ¼$sfðs;gÞ, with $equal to the gradient operator with respect to the variables . The evolution of the elastic left Cauchy–Green tensor (i.e. flow rule) is defined in the ALE setting as [14] obe otjvþc$xbelbebelT¼2_ ccmðs;gÞbe;ð6Þ where l¼$xvis the velocity gradient tensor and _ cc is the plastic multiplier. Equation (6) is complemented with the Kuhn-Tucker conditions _ cc 0;fðs;gÞ0 and _ ccfðs;gÞ¼0:ð7Þ The evolution of the relative density gis given by the mass conservation principle, which in ALE formulation reads og otjvþc$xgþg$xv¼0:ð8Þ 2.3 Numerical time-integration: a fractional-step method In the ALE approach, the two fundamental unknowns in every time-step ½nt;nþ1t, with nþ1t¼ntþDt, are the increment of mesh displacements, nþ1DuUðvÞ¼nþ1xðvÞnxðvÞð9Þ and the increment of particle displacements, nþ1DuðnXðvÞÞ ¼ nþ1xðnXðvÞÞ  nxðnXðvÞÞ ð10Þ which is referred to the particles nXassociated to grid points vat the beginning of the time-step, nX¼Xðv;ntÞ. The mesh displacements are obtained from an ALE remeshing algorithm. The particle displacements are found via the incremental solution of the equilibrium equation. The difference between the increments of particle and mesh displacements is the so-called increment of convective displacements, nþ1DuconvðvÞ¼nþ1DuðvÞnþ1DuUðvÞ;ð11Þ which represents the relative motion between particles and grid points during the time-step. The quantities to integrate are only beand g, see remark 1 below. The numerical time-integration is done by means of a fractional-step method over Equations (6) and (8). Every time–step is divided into two phases: the Lagrangian phase and the convection phase. During the Lagrangian phase convection is neglected, c¼0. In that situation, Eq. (6) simplifies to the standard spatial expression of the flow rule Lvbe¼2_ ccmðs;gÞbe;ð12Þ where Lvbe¼obe otjxlbebelTis the Lie derivative with respect to the particle velocity; and Eq. (8) simplifies to gðX;tÞ¼g0ðXÞ detðFÞ;ð13Þ the material expression of the Lagrangian relative density evolution. The initial conditions for equations (12) and (13) are the state at time nt(nx,nbeand ng). The problem is strain-driven, thus, the Lagrangian values Lbeand Lgare obtained from nþ1Du. Recall that nþ1Duis not required to be in equilibrium. After the Lagrangian phase, an ALE remeshing algorithm is employed to compute the increment of mesh displacements nþ1DuU. The algorithm is defined such as it reduces the element distortion of a pure Lagrangian approach. During the convection phase, the convective term of the evolution equations is taken into account. The convective velocity is assumed constant in the time-step and equal to c¼nþ1Duconv=Dt. The final (convected) values nþ1beand nþ1gare computed integrating the transport equations obe otjvþc$xbe¼0and og otjvþc$xg¼0;ð14Þ with the initial conditions Lbeand Lg. In this work, only the Lagrangian phase is performed during iterations for equilibrium within each load increment. The remeshing and the convection are computed with the converged results at the end of each time-step. This results into a reduced computational cost overhead with respect to the pure Lagrangian approach. The equilibrium perturbation due to the remeshing and the convective time-integration of the converged results is taken into account in the subsequent load increment. Therefore, there are no accumulative errors in the results. On the other hand, as the iterative scheme is pure Lagrangian, the consistent operator is not influenced by the remeshing technique nor by the numerical time-integration scheme used in the convective phase. The overall scheme is summarized in table 1. The main characteristics of the numerical time-integration of each phase are presented in the following two subsections. Remark 1 In a general setting, all quantities related with the particles should be convected. This includes the material parameters and the value of detðFÞ, in addition to be and g(and the plastic internal variables if they exist). However, in the usual case that the material surfaces are tracked by the ALE remeshing algorithm and homogeneous materials are considered (as is the case in all the examples presented in this work), only the integration of beand gis necessary. The value of detðFÞcan be computed from the material expression of the relative density evolution, Eq. (13). 222 2.3.1 Lagrangian phase The results of the Lagrangian phase, Lbeand Lg, are computed from nbe,ngand the incremental particle deformation gradient, Lf¼Indim þr nxnþ1Du;ð15Þ which relates the particle deformation gradients at time nt, nF, and at the end of the Lagrangian phase, LF, through the relationship Lf¼LFðnFÞ1.Indim denotes a identity matrix of order ndim. The relative density, g, is integrated exactly (in the sense that no numerical time-integration scheme is used) because of the material expression of the mass conservation principle, Eq. (13), which leads to Lg¼g0 detðLFÞ¼ ng detðLfÞ:ð16Þ The value of Lbeis obtained by means of the standard elastic predictor-plastic corrector split strategy applied to Eqs. (7) and (12). Remarkably, the dependence of the constitutive equations on the density does not modify the algorithm. The value of Lgis given by Eq. (16), and therefore it plays the role of a fixed parameter. The result of the elastic predictor step is the so-called trial state. It is defined by trbe¼LfnbeLfT:ð17Þ If the trial state is admissible, fðsðtrbeÞ;LgÞ0, the value of Lbeis set equal to the trial state. If it is not, a plastic corrector step is computed. Recall that the trial state, Eq. (17), does not depend on Lg, a key point in the linearization of the algorithm presented in Sect. 3. The plastic correction step requires the approximation of the flow rule, Eq. (12). A standard approximation consists in the use of the exponential map and the backward Euler integration scheme [3]. Under the previous isotropy assumptions, this approach leads to a nonlinear system of equations with the same structure as that of infinitesimal elastoplasticity. In order to obtain this nonlinear system of equations, three vectors of Rndim are defined: tree,Leeand L ss. The components of treeand Leeare functions of the eigenvalues of the tensors trbeand Lbe, respectively, and the components of L ss are the eigenvalues of Ls: *be¼X ndim i¼1expð½eeiÞ2ni tr ni tr and Ls¼X ndim i¼1 ½L ssini tr ni tr ; ð18Þ where the superscript *refers to tr and L, fni trgi¼1;...;ndim are the eigenvectors of the three tensors, and f½*igi¼1;...;ndim the components of the three vectors. The eigenvectors are the same for the three tensors because of the isotropy assumptions. Thus, they are fully specified by the trial state. After some manipulations, the following nonlinear system of equations is found: LeeþDcmss  ssðLeeÞ;Lg  ¼tree f ssðLeeÞ;Lg  ¼0;ð19Þ where Dc¼Dt_ cc is the incremental plastic multiplier, mss is the flow vector in the principal direction space (that is, ms¼Pndim i¼1½mssini tr ni tr), the isotropic functions mss and fare expressed as functions of  ss, and  ss is given by the hyperelastic relationship  ss ¼d WWe dee;ð20Þ with  WWeðeeÞdefined so that Eq. (20) is equivalent to Eq. (4). Equations (19) are complemented with the restriction Dc0. This problem is typically solved with the Newton– Raphson method. Remarkably, as Lgis fixed, the globally convergent extensions of the Newton–Raphson method presented [22] for density-independent finite strain models also apply to density-dependent plastic problems. Once Eqs. (19) are solved, the Lagrangian phase is completed. 2.3.2 Convective phase During the convective phase, the transport equations (14) are integrated numerically with a Godunov-like scheme [13, 24]. The application of this scheme to quasi-static problems is compared with other alternatives in [13]. This scheme circumvents the computation of gradients with respect to the spatial coordinates present in the transport equations, the main source of trouble when discretized problems are considered. The only difference between the present application, with hyperelastic-plastic models, and that presented in [13, 24], with hypoelastic-plastic models, is that here it is applied to the convection of the elastic left Cauchy–Green tensor instead of the Cauchy stresses. The key points of the implementation are presented in the following. The transport equations (14) contains seven equivalent scalar equations, one for each component of beand one for g. Therefore, the scheme is devised for a general scalar convective equation ou otjvþc$xu¼0;ð21Þ Table 1. The overall ALE scheme For every time-step [ n t, nþ1 t]: Material phase Neglect convective terms Advance the solution iteratively in an updated Lagrangian fashion: compute the increment of particle displacements n+1 Du and quantities L b e and L g(superscript L denotes Lagrangian) Remeshing Compute the increment of mesh displacements n+1 Du F and the increment of convective displacements n+1 Du conv by means of a remeshing algorithm that reduces element distortion. Compute the convective velocity c= n+1 Du conv /Dt Convection phase Account for convective terms Use the Godunov-type technique to convect quantities L b e and L ginto n+1 b e and n+1 g Compute stresses n+1 sand n+1 r 223 with uðt;xÞrepresenting the different components of be and g. The initial condition for Eq. (21) is the Lagrangian value Lu. The basic idea of the scheme is to divide every finite element in various zones, each of them corresponding to the influence domain of a Gauss point. In two-dimensional problems, as those presented in this work, each zone has an area Aand nedg edges. Then, at each Gauss point, the following explicit update equation is applied nþ1u¼LuDt 2AX nedg iedg¼1 Iiedg Luc iedg Lu  1signðIiedg Þ; ð22Þ where Luc iedg is the value of the variable uin the contiguous Gauss point across edge iedg and Iiedg is the flux of convective velocity cacross edge iedg. This scheme leads to a very simple algorithm. Moreover, as the space discretization and the convective velocity are the same for the different scalar transport equations, the major part of the computations are common to all of them. Once Eq. (22) is applied to the six components of beand gthe convective phase is completed: the values of nþ1be and nþ1gare determined. The stress values nþ1sand nþ1r are computed with Eqs. (4) and (5). 3 Consistent tangent operator The consistent tangent moduli are needed to solve the equilibrium equation with quadratic convergence [20]. In this work, as the convective phase is not included within the iterative process, they are the same needed in standard Lagrangian approaches. However, due to the densitydependence of the plastic equations, the well-known consistent tangent moduli for density-independent models for multiplicative finite strain problems [3] results incomplete. In the following, the complete expression for densitydependent plastic models is presented. The consistent tangent moduli, c, are the linearization of the Kirchhoff stresses obtained from the Lagrangian phase, Ls, with respect to the gradient of the incremental particle displacements, rnxnþ1Du[3, 25]. They are found by applying the chain rule to equation (18)2, which leads to c¼X ndim i¼1X ndim j¼1 ½aijni tr ni tr nj tr nj tr þ2X ndim i¼1 ½L ssi^ cci tr ; ð23Þ with the first term corresponding to the linearization of L ss and the second one to the linearization of fni trgi¼1;...;ndim , and where the tensors f^ cci trgi¼1;...;ndim depend on fni trgi¼1;...;ndim and tree[3], and ais a matrix of order ndim defined as a¼dL ss dtree:ð24Þ Equation (23) has the same expression for densitydependent and density-indepedent plastic models because the density does not affect the application of the chain rule nor the trial state. In fact, the influence of the density is restricted to the values of L ss and ain plastic steps, i.e. when the trial state is not admissible. In elastic steps, the matrix ais the Hessian of  WWe, a¼d2 WWe dee2ee¼tree :ð25Þ The value of L ss is obtained, in both elastic and plastic steps, directly from the hyperelastic relationship. Therefore, only the expression of afor plastic steps needs to be determined. In order to do that it is useful to rephrase the dependence of Lgon Lf, Eq. (16), as Lg¼n^ gg exp trðtreeÞ;ð26Þ with n^ gg ¼ngffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi detðnbeÞ pa known value from the previous time step and where trð*Þmeans trace of (that is, Pndim i¼1½*i). Equation (26) is found by applying the determinant function to both sides of Eq. (17) and substituting Eq. (18)1into it. A more convenient expression of ais found from its definition, Eq. (24), and the application of the chain rule: a¼dL ss dLee dLee dtree¼d2 WWe dee2ee¼Lee dLee dtree:ð27Þ This expression shows that ais determined once the total influence of treeon Lee;dLee dtree, is found. This influence is given by the the nonlinear system of Eq. (19) and the relationship between Lgand tree, Eq. (26). Thus, it can be computed by linearizing the nonlinear system of equations LeeþDcm ss  ssðLeeÞ;Lg  ¼tree f ssðLeeÞ;Lg  ¼0 Lg¼n^ gg exp trðtreeÞ: ð28Þ Once dLee=dtreeis determined, the expression of ais found substituting it into Eq. (27). The result can be rearranged as a¼agindep þag;ð29Þ with agindep equal to the standard consistent moduli for density-independent plastic models, agindep ¼~ GG  ~ GGm ss r ssfT~ GG r ssfT~ GGmss ;ð30Þ where ~ GG ¼d2 WWe dee21 þDcr ssm ss1 ;ð31Þ and ag¼gDc~ GG  ~ GGm ss r ssfT~ GG r ssfT~ GGm ss om ss og11;ndim þg of og ~ GGm ss r ssfT~ GGm ss 11;ndim ;ð32Þ a term that takes into account the influence of the density. All quantities involved in the computation of a, Eqs. (30–32), are evaluated at the end of the Lagrangian phase. 224 In the density-independent case, symmetric tangent moduli are obtained for associative material models. On the contrary, unsymmetric moduli are found with all density-dependent material models because agis, in general, unsymmetric. For this reason, in density-dependent plasticity, unsymmetric linear solvers have to be used in order to keep the characteristic quadratic convergence of the Newton–Raphson method. On the other hand, it is important to remark that the expression of agcan be computed with just a few more matrix–vector products than the standard agindep, compare Eqs. (32) and (30), and, as expected, the additional information om ss ogand of og:ð33Þ 3.1 Numerical differentiation A key point in the computation of consistent tangent moduli in density-independent plasticity is the computation of flow vector and hardening law derivatives. Reference [18] presents an analysis of different numerical differentiation schemes to compute these derivatives and references [23, 26] show the application to different nontrivial elastoplastic models and time-integration rules. The main conclusion of these works is that simple first and second order difference schemes do not disturb the quadratic convergence of the Newton–Raphson method provided that the stepsize is fixed in a relative way. Moreover, in that references it is shown that numerical differentiation is an efficient alternative to analytical derivatives even if are readily available. These schemes are directly applicable to densitydependent plastic models. In this case, it can also be applied to the computation of the derivatives of flow vector and the yield function with respect to the relative density, see Eq. (33), a situation of special interest when realistic constitutive laws are considered. In this work, a first order difference scheme with a relative stepsize equal to 105has been used to approximate these derivatives. 4 Numerical simulations In this section, the proposed approach is applied to the simulation of several powder compaction problems. In the first two cases the attention is focused in the comparison of the present numerical results with numerical and experimental ones reported in the literature. The first one is solved with a pure Lagrangian approach, and the second one using the ALE formulation presented in Sect. 2. Specific analysis related with the application of the ALE formulation, the definition of the globally convergent schemes for the plastic corrector step, and the computation of the consistent tangent moduli via numerical differentiation schemes can be found in references [14, 18, 22, 23]. The powder material is modelled with the Hencky’s hyperelastic law, which leads to a linear relationship between  ss and ee, and associate plasticity with the following elliptic yield function [1, 27] fellipðs;gÞ¼2J2ðsÞþa1ðgÞI1ðsÞ 3  2 2 3a2ðgÞðryÞ2; ð34Þ with I1ðsÞequal to the first invariant of s;J2ðsÞequal to the second invariant of the deviatoric part of s, and the density–dependent parameters a1ðgÞ¼ 1g2 2þg2  n1 g<1 0g1 (and a2ðgÞ¼ 0:02g0 10:98g0  n2 gg0 g0:98g0 10:98g0  n2 g>g0 8 > < > : ð35Þ The dependence of a1and a2on gfor the material parameters presented in Table 2 is depicted in Fig. 2. The trace of the yield function on the meridian plane psqs, with ps¼I1ðsÞ 3and qs¼ffiffiffiffiffiffiffiffiffiffiffiffi 3J2ðsÞ p, for different relative densities are depicted in Fig. 3. Note that fellip becomes the von Mises yield function for g1. 4.1 Compaction of a plain bush component The first example is the uniaxial compaction of a plain bush component. Experimental data [28] and numerical results [2, 28, 29] are available for this example. Both are used here for comparative purposes. The test is performed with the powder material parameters presented in Table 2. These material parameters are calibrated in [29] by comparing the numerical results with those presented in [28]. The component is modelled by an axisymmetric representation as illustrated in Fig. 4 [2]. A 2D structured mesh of 200 bilinear elements is used. The die wall friction is simulated with a Coulomb friction coefficient l¼0:15 Table 2. Material parameters E50 000. [MPa] m0.37 r y 12. [MPa] g 0 0.41 g 1 0.5 g 2 2.2 Fig. 2. Dependence of parameters a1ðgÞand a2ðgÞon the relative density, g 225 [28, 29] acting in the inner and outer walls of the sample, segments BC and DA of Fig. 4. In [29] the relative radial movement of top and bottom surfaces with respect to the punches, segments AB and CD, is allowed, with a Coulomb friction coefficient equal to that of the lateral surfaces. Here, following [2], the radial displacement of the top and the bottom of the sample is restrained. The vertical displacement of segment AB is also set equal to zero, and a vertical displacement of 11.5 mm is imposed to segment CD, simulating the top punch movement. The relative density profile at radius equal to 10.5 mm for a top punch displacement of 10 mm is shown in Fig. 5. The numerical results of [28, 29] and the experimental data of [28] are included in the same figure. The results of the present numerical simulation are in agreement with all of them. However, as expected because of the similarities of the material model, the best agreement is found with the results of [29]. The lower relative density at the top part of the sample found by [29], with respect to the present one, and the higher values in the lower part are directly related with the different treatment of the friction effects of top and bottom punches. The evolution of the vertical reaction of the punch with respect to its vertical displacement is depicted in Fig. 6. The numerical results of [2] and the experimental data of [28] are also included. In reference [2] the results obtained with different formulations are shown. Those depicted in Fig. 6 are the most similar to the experimental data. The agreement between the results of the present simulation and the experimental ones is good. Finally, the evolution of the relative density distribution over the sample is depicted in Fig. 7. Results for different vertical top punch displacements (from 4 to 11 mm) are shown. Two different relative density scales are used (one in each row) in order to show better the non-homogeneous distribution of the density. The shape of the relative density distributions are in general agreement with those Fig. 3. Trace of the elliptic yield function on the meridian plane qspsfor different relative densities, g Fig. 4. Plain bush component. Problem definition (after Lewis and Khoei 1998) and computational mesh Fig. 5. Plain bush component. Relative density profile at radius equal to 10.5 mm Fig. 6. Plain bush component. Relationship between top punch vertical reaction and its vertical displacement 226 presented by [2], although the quantitative values are a little bit higher. In order to check the present results, the evolution of the sample mass during the simulation has been computed. A constant value equal to the initial one has been found, up to six significant digits. Therefore, the mass conservation principle is verified. This implies, for instance, a mean value of the relative density equal to 2g0for a height reduction of 50% (i.e. gmean ¼0:82 for a top displacement of 10 mm, see Figs. 6 and 7). 4.2 Compaction of a rotational flanged component The second example correspond to the compaction of a flanged component which is modelled by an axisymmetric representation, as illustrated in Fig. 8. The example is used in [2] to illustrate the applicability of a dynamic approach with a hypoelastoplastic model. The same example is used in [6] to show the utility of an h-adaptive remeshing technique to reduce mesh distortion. Here the hyperelastoplastic model and the ALE formulation presented previously are used. Some experimental results [30] are available. The present results illustrate that the proposed approach allows to simulate highly demanding powder compaction processes without mesh distortion and spurious oscillations in the results. Moreover, it is shown that the mass conservation principle is verified with a low relative error. Three different compaction tests are simulated [2]: 1) a vertical movement of the top punch (6.06 mm); 2) a vertical movement of the bottom punch (5.10 mm); and 3) a simultaneous movement of both punches (6.06 mm the top punch and 7.70 mm the bottom punch). The same structured mesh of 170 eight-noded elements with reduced integration (four Gauss points per element) is used in the Fig. 7. Plain bush component. Relative density distribution for different top punch movements. Note that two different scales are used 227 three tests. The die wall friction is simulated with a Coulomb friction coefficient l¼0:08 acting in the segments BC, CD, DE and FA, see Fig. 8, and the radial displacement at the punches is restrained [2]. The analysis is performed with the powder material of Table 2, calibrated in [1] for the compaction of the plain bush component. The relationships between dimensionless vertical loads and punch displacements obtained in this work are compared with those of [2] in Fig. 9. The agreement between the two sets of curves is evident. However, the reference load is different: 392 kN for the present results and 1550 kN for the results of [2] (in both cases the final reaction of the top punch in the doublepunch compaction test). This difference can be related with the different modelling approach (dynamic versus static) and, especially, with the great difference between the elastic moduli used in both cases (40 MPa in [2], three orders of magnitude lower than the one used here). In the following, a detailed analysis of the relative density distribution obtained in each of the three tests is presented. In the first test, top punch compaction, the mesh region ABCG is Eulerian, and equal height elements are prescribed in the mesh region GDEF. The variation of the mass of the sample during the simulation is depicted in Fig. 10. A final loss of 0.75% is found. This small variation corresponds to the truncation error in the discretization of the convective term of the ALE formulation (both temporal and spatial discretization). The evolution of the relative density distribution is summarized in Fig. 11(a–d). The compaction process leads to a clearly non-homogenous density distribution. As expected, higher values are found in the outer region of the sample and lower ones close to the bottom surface. A smooth transition from higher to lower densities is found. A dense zone is detected in the corner region, just over the point C, during all the process. Recall that there is not mesh distortion because, although the flux of mass is important, the mesh does not follow the material particles in the ALE formulation. The final relative density profile at 1.88 mm from line GD is depicted in Fig. 11(e). The present results are in general agreement with those presented in [2]. Two zones with a quasi-uniform relative density are found in both cases. However, the density profile obtained in [2] presents a big oscillation between these two zones, and, on the contrary, the transition obtained in this work is smooth. The second test consists in a bottom punch compaction. In this case the mesh region GDEF is Eulerian and equal height elements are prescribed in region ABCG. The variation of the mass of the sample during the simulation is depicted in Fig. 10. The mass gain is less than 1% at the end of the simulation. The error is a bit larger than for the top punch test. This indicates that the convective effects are more important in this case. Fig. 8. Flanged component. Problem definition (after Lewis and Khoei 1998) and computational mesh Fig. 9. Flanged component. Relationships between the vertical reactions of the punches and their vertical movements: top reaction for top punch compaction, bottom reaction for bottom punch compaction, and top and bottom reactions for doublepunch compaction Fig. 10. Relative mass variation during the three compaction processes of the flanged component. The load levels are referred to the punch displacements imposed at the end of each test 228