Regularization of shear banding and prediction of size effects in manufacturing operations: A micromorphic plasticity explicit scheme

Good quality manufacturing operation simulations are essential to obtain reliable numerical predictions of the processes. In many cases, it is possible to observe that the deformation localizes in narrow areas, and since the primary deformation mode is under shear, these areas are called shear bands. In classical continuum mechanics models, the deformation localization may lead to spurious mesh dependency if the material locally experiences thermal or plastic strain softening. One option to regularize such a non-physical behavior is to resort to non-local continuum mechanics theories. This paper adopts a scalar micromorphic approach, which includes a characteristic length scale in the constitutive framework to enforce the plastic strain gradient theory to regularize the solution. Since many manufacturing process simulations are often assessed through finite element methods with an explicit solver to facilitate convergence, we present an original model formulation and procedure for the implementation of the micromorphic continuum in an explicit finite element code. The approach is illustrated in the case of the VPS explicit solver from ESI GROUP. According to the original formulation, we propose an easy way to implement a scalar micromorphic approach by taking advantage of an analogy with the thermal balance equation. The numerical implementation is verified against the analytical solution of a semi-infinite glide problem. Finally, the correctness of the method is addressed by successfully predicting size effects both in a cutting and a bending tests.


Introduction
It is well known that the classical Cauchy continuum description is not sufficient to predict the different responses of the medium when either stresses or strains localize. Although experimental evidence strongly emphasized the existence of size-dependent behaviors, where smaller is the size, stronger is the response, the classical continuum mechanics models do not possess a characteristic length scale that allows the prediction of said size-effects. The description of the classical continuum mechanics is, in fact, of a local nature, meaning that the configuration of the medium at any location is solely determined by the properties characterizing the continuum at that specific location, and the distribution of the said properties in the neighborhood of this location does not influence the local properties. Several descriptions of the continuum have been proposed in the literature as alternatives to the classical continuum mechanics, in the attempt of including a gradient-related response of the medium, and these theories are referred to as non-local or higher-order theories. In general, any continuum mechanics model, different from the classical model, belongs to the family of the generalized continuum mechanics. A common feature that is shared by all of them is the appearance of a characteristic length scale in the constitutive framework, which naturally arises when the internal power and the constitutive material model are explicitly defined [25].
Another limitation of the classical theory can be identified when tackling strain localization problems: if a specific form of the constitutive material behavior is chosen, the static boundary value problem loses its ellipticity and assumes a hyperbolic character. This problem characterizes, but is not limited to (same behavior could be found for non-associated plastic flow, e.g., [29,37,43]), the condition in which the material tangent experiences a local negative slope. The change in the form of a boundary value problem causes the solution not to be uniquely determined anymore. The negativeness of the material tangent could be induced by several physical phenomena, and the one which most concerns our research is the thermal softening induced in metallic materials severely deformed. Such behavior can be experienced by the continuum when high temperatures are locally produced by the plastic deformation and subsequently retained due to the the combination between low thermal conductivity of the materials and high strain rates. At relatively high temperatures, the material reduces its yield strength, subsequently experiencing a softening of a thermal nature. From the analytical point of view the boundary value problem is not uniquely defined, and from the numerical point of view (if the problem ought to be discretized through Finite Element Method, for instance), the solution appears to be spuriously mesh-dependent.
These problems are severely relevant to manufacturing process simulations because the material is heavily deformed in a short amount of time, thus inducing deformation localization and thermal softening. Moreover, it was already demonstrated that a strong size effect characterizes some manufacturing processes and that the classical continuum mechanics was no longer adequate to predict these behaviors [7,16,24,31,44].
The present contribution aims to assess these two main problems (size-dependency and mesh-dependency) in metal manufacturing operations. The main objective is to propose a method that is capable of solving both problems. Resorting to a generalized continuum mechanics theory could compensate for the limitations of the classical continuum mechanics both in terms of size-effects and spurious mesh-dependency. Several contributions can be found in the literature focusing on the description and comparison among many different non-local theories, some of them focusing on the so-called strain gradient theories [15,19], or the gradient of internal variables [20], or elastic-gradient theories [3], or specifically exploiting the topic of using a generalized continuum mechanics theory for manufacturing operations [34].
In recent years, different generalized continuum mechanics theories have been used to simulate manufacturing operations. A flat punch molding process was simulated by Guha and co-workers using the plastic strain gradient theory [16]. A similar theoretical framework was used to simulate steady-state rolling processes in [30]. Other relevant contributions in the application of non-local theories for manufacturing operation simulations can be found in the investigations of Liu et al. [22][23][24] and Asad et al. [2], who both reproduced orthogonal cutting simulations using strain gradient effects. More recently, Diamantopoulou and coworkers [8] used a non-local continuum mechanics theory enhanced with the gradient of a scalar damage variable to simulate metal forming.
Besides the strain gradient theories, other generalized continuum mechanics theory could solve the aforementioned problems. The micromorphic approaches, for instance, involve the gradient of a tensor of any rank, supposed to perform the targeted strain gradient operation [12], and they can also be used to overcome the limitations of the classical continuum mechanics. The micromorphic approach introduces additional degrees of freedom in the problem, and, depending on the type of theory that is required, the computational cost might dramatically increase. For instance, a full-order micro-curl model, as the one proposed by Cordero et al. [5], requires at least 5 additional degrees of freedom in two-dimensional settings. In contrast, the micromorphic approach involving a scalar micromorphic variable, so-called reduced-order micromorphic model, includes only one additional degree of freedom.
In this context, our contribution aims at investigating the size-effect predictions and regularization properties of a time-dependent strain gradient theory that is implemented through a scalar micromorphic framework using an explicit formulation, in which a viscoplastic micromorphic-related variable is included, but no micromorphic inertia is present. The main novelty of the proposed method lies in the easiness of the implementation of the theory in an already-well-structured finite element solver. It represents an alternative to existing implementations of such micromorphic models like the one proposed by Saanouni and Hamed [36]. The framework that we will present can, in fact, simply be solved through a common thermal-field solver, and such crucial aspect will be properly addressed in the present paper.
The layout of the manuscript is as follows. The formulation of the analytical model is provided in Section 2 in which both the kinematics and the energetic aspects of the theory are presented, alongside its thermodynamic description, so that the recoverable and dissipative contributions are explicitly stated as such. The section concludes with the pivotal analogy between the thermal and the micromorphic balance equations, which further simplifies any possible implementation of the theory in a Finite Element software. In Section 2.7 the discretization of the equations and the implementation of the method in the explicit finite element software VPS/Pam-Crash ® from ESI [10] will be described. Section 3 will be used to present a simple analytical solution that will be useful to verify the implementation of the model in a finite element framework. Finally, in Section 4 the numerical method will be used to simulate two manufacturing operations in which significant strain gradient effects are expected to take place, namely the shear/trimming operation and the bending test. The mesh-dependency will be analyzed, along with the size-effect in terms of cumulative plastic strain distribution. Conclusions follow in Section 5.
Notations In this work, the following notations are used. The first, second and fourth-order tensors will be indicated by a bar, tilde or double tilde, respectively, underneath the tensor: a, a ∼ and a ≈ . The single or double dot above of a degree of freedom indicates the first or second time derivative respectively. Single and double contractions are indicated by single dot or double dots between the tensors and they operate on the inner indices of the tensors: where Einstein index summation convention applies. The single and double tensor products operate as:

Theoretical formulation and finite element implementation
The micromorphic scheme has been proven to be a straightforward and relatively simple procedure to introduce additional degrees of freedom to the continuum in order to achieve non-local regularization effects [12,13], and it has been used already in several other contributions [1,6,8,26,32,35,36]. Among the cited works, the only ones to adapt and implement the micromorphic approach for an explicit time-dependent problems can be found in [6,8,35,36]. The aforementioned authors presented a time-dependent framework, in which the governing equations for the micromorphic variables include a second-order time derivative of the micromorphic variables. Additional coefficients associated with this term were included to characterize the inertia of the micromorphic variables, a role that is usually assigned to the density for the governing equations of displacement fields. Furthermore, Davaze and co-workers [6] included some dissipation terms associated with the first-order time derivative of the micromorphic variable in governing equation so as to avoid any oscillation of the solution caused by the form of the partial differential equation (specifically induced by the presence of a second-order time derivative term). However, the authors of this research used the theory to achieve mesh-regularization for fracture growth simulations in metals. Exploring the extent of such an approach for manufacturing operation simulations was not their target.
In the present work, we make use of a scalar micromorphic approach to govern the strain gradient effect and to restore mesh independence. The classical continuum mechanics model is enhanced with one additional degree of freedom. The governing equations for such an additional variable will be directly derived from the definition of the internal power. The micromorphic approach will be used to control the distribution of the cumulative plastic strain. Therefore, the additional degree of freedom will be enforced to follow this quantity through a penalty term.
In this section, the kinematics of the theory will first be provided, from which the balance equations can be derived, the definition of the Helmholtz free energy and of the Clausius Duhem inequality will follow. Finally, the section will conclude with the analogy between the micromorphic-balance equation and the thermal field equation.

Kinematics and balance equations
The kinematics of the model follows the one commonly used in the classical continuum mechanics. The secondorder strain tensor is defined as: with u being the displacement vector and denotes the gradient of a vector. Furthermore, the total strain tensor is additively decomposed into an elastic part e ∼ and a plastic part p ∼ as follows: By indicating the velocity v as u , we can define the strain rate as: The variables which are supposed to carry the targeted strain gradient effects are selected among the available state variables which can be tensors of any rank. Here, we consider a scalar variable. Two types of degrees of freedom (DOFs) are applied to the material point: the classical displacement vector u and the additional scalar micromorphic variable p associated with the cumulative plastic strain p through the penalty term H to be defined later. Therefore, every node is endowed with 3 displacement and 1 micromorphic variable: Based on the definition of the strain and of the micromorphic variable, we are allowed to write the internal and kinetic power densities of the body as dependent on the strain, the micromorphic variable and its gradient 1 : where is the mass density and ü is the acceleration vector. Here, a and b are generalized stresses associated with the micromorphic variable and its gradient, respectively. In this formulation, the densities of power generated by external forces and contact forces can be written as: with f e being the density of body force, a e and b e are the generalized body stresses associated to p and its gradient. f c and a c are the classical traction and the micromorphic traction. The contact power density defined in Eq. 8 clearly states that the gradient of the micromorphic variable is not linked to any boundary effect. The global power balance law can be written as: which, through Eqs. 5, 6, 7 and 8, transforms into: p (e) =f e ⋅u + a eṗ + b e ⋅ ṗ , Based on the principle of virtual power, the equilibrium equations are obtained as: which are bounded by the following Neumann boundary conditions: where n is the outer normal to the surface closing the domain .

Helmholtz free energy potential
The constitutive model of the medium characterizing the shape of both the classical and the generalized stresses is provided via the definition of their associated potential. The free energy density function is assumed to depend on the following state variables: namely, the elastic strain, the cumulative plastic strain, the micromorphic variable, and its gradient. The chosen potential has the form: where C ≈ is the elastic fourth-order stiffness tensor, p is the plastic contribution to the Helmholtz free energy (in case of hardening/softening it accounts for the expansion/shrinking of the yield surface in the stress space), and is the additional micromorphic contribution. A simple quadratic potential is adopted for the latter: where A ∼ is the higher-order modulus. For an isotropic material, the elastic stiffness tensor C ≈ and the tensor of higherorder moduli reduce to the following forms: 1 There is a possibility here to explicitly define the kinetic and damping energy of the continuum as function of the micromorphic variable as well. Such type of descriptions have already been proposed by other researchers [6,28,36]. In the present work, however, we will include instead a viscous contribution of the micromorphic variable in the constitutive model of the continuum. with and as the classical Lamé parameters and A the new higher order modulus. The following linear isotropic plastic behavior is assigned to the material: where H p is the hardening modulus. Nonlinear hardening laws are possible but not considered here for simplicity.

Clausius-duhem inequality
The Clausius-Duhem inequality will be used to ensure the thermodynamic consistency of the model and to define the recoverable and dissipative parts of the mechanical contributions. The local form of the second law of thermodynamic for an iso-thermal transformation can be expressed for a continuum body as: Expanding the time derivative of Helmholtz free potential with respect to the variables on which it depends, and by retrieving the additive elasto-plastic decomposition of the strain rates, the Clausius-Duhem inequality reads: At this stage, the choice on the elastic part of the strain to be energetically recoverable can be made. It implies that the terms multiplying ̇ ∼ e must vanish so as to ensure that the elastic strain does not contribute in entropy production, leading to: The distinction between recoverable and dissipative parts of the generalized stress terms must also be drawn. For the gradient of the micromorphic variable, we assume that it is fully recoverable, therefore: This means that the gradient of plastic strain solely contributes to the free energy potential. In the case of metals, this can be justified by the fact that the plastic strain gradient contains contributions of the dislocation density tensor which is known to be associated with energy storage [14,18]. Regarding now the dissipation produced by the variation of the micromorphic variable, its positiveness can be ensured, as originally suggested by Gurtin [12,17], by imposing that the generalized stress a possesses a recoverable part and a dissipative part that depends on ṗ itself: where C is a parameter related to viscous micromorphic effects. Lastly, for the plastic part of the Helmholtz free energy: where R is a thermodynamic force associated to variation of the cumulative plastic strain. The residual dissipation rate can now be written as: The positiveness of the new parameters A and C then ensures the positive definiteness of the micromorphic contributions in the free energy density and in the dissipation rate.

Partial differential equation governing the micromorphic variable and enhanced hardening law
By considering the explicit definition of the Helmholtz free energy potential given in Eq. 17, the generalized stresses read: The previous Eq. 28 indicates that the micromorphic variable p and cumulative plastic strain p are related to each other through the penalty term H . In order for the micromorphic variable to closely match the value of the cumulative plastic strain, it is necessary to ensure that the value of H is relatively large. At this stage, it is possible to re-write the additional partial differential equation governing the micromorphic distribution by plugging the selected constitutive behavior into it. In absence of higher-order body forces ( a e and b e ), Eq. 12 can be written as: where ∇ 2 indicates the Laplacian differential operator. The previous equation represents the only additional equation that must be solved combined with the ones governing the displacement fields. Previous researchers already explored the potential of the micromorphic theory in rate-dependent analysis under explicit integration schemes using a modified version of Eq. 30. For instance, Saanouni and Hamed proposed a theory in which the second-order time derivative (acceleration) of p takes the place of the first-order time derivative in Eqs. 30 [8,36]. Therefore, in analogy with the PDE governing the displacement fields, a form of inertia was associated to the micromorphic variable, whereas, in case of the present investigation, a viscous term associated to the micromorphic variable is considered. The PDE governing the micromorphic field can be rewritten as: where l ch is the characteristic length scale endowing the theory with the spatial regularization property, and ch is a characteristic time. To fully solve Eq. 31, it must be coupled with a constitutive model for the plastic behavior of the medium. Starting from the yield function: where eq is the von Mises equivalent stress measure and 0 is the initial yield stress. Assuming associated plasticity and the normality rule to hold, the rate of the plastic strain can be written as: and the dissipation in Eq. 27 takes the form: and in case of plastic loading: From the specific form of the plastic part of the Helmholtz free energy and from Eq. 26, one can infer the stress that is thermodynamically associated with the cumulative plastic strain: This represents the hardening law enhanced by a new micromorphic contribution. This shows the coupling arising in the theory between plasticity and the micromorphic variable. After substituting the Eq. 31, the following alternative expression of the enhanced hardening law is obtained: The linear hardening/softening contribution (depending on the sign of H p ) to the yield stress is enhanced by the Laplacian of the micromorphic variable, a usual term in regularization methods, but also by an additional viscous term whose magnitude is controlled by the value of parameter C .

Micromorphic-thermal analogy
The comparison between the scalar micromorphic model described in the previous section and the classical thermomechanical theory is here outlined. The development of the latter theory will not be fully reported, but we will make use of the main governing equations of the thermal field to draw the comparison with the micromorphic theory previously developed. On the one hand, the additional variable in the present theory, p , ought to be solved through the PDE Eq. 30, whereas, on the other hand, the additional degree of freedom of the classical thermomechanical theory, that is temperature T, must be solved through a different PDE, and here the two equations are reported (where the Fourier conduction law is assumed to be valid for the heat flux) where C is the specific heat capacity of the material, r is a source term and k is the thermal conductivity of the material, that we assumed to be independent from temperature. Although the two equations are used to govern completely different physical fields, a straightforward parallelism among them can be identified. In Table 1, a comparison between different aspects of the two theories is reported. The analogy between these two theories inspired the idea of adapting an already implemented numerical resolution scheme (meant to be used for the thermal field) for the micromorphic variable. The main objective of the present investigation is, in fact, the analysis of the feasibility of such idea. The main advantage of the proposed method is that the micromorphic theory can be easily implemented in an explicit resolution scheme, while requiring very limited access and marginal effort in modifying the original code. This aspect obviously makes the implementation of this theory more attractive than other methodologies which would require high level of accessibility to the main solver, since both new element and material definitions would need to be developed. Such an analogy has been used in the past for coupling chemical diffusion and mechanics in the implicit version of the numerical solver ABAQUS [9]. The analogy has also been recognized and used to implement gradient plasticity and gradient damage models, again, in the implicit version of ABAQUS [40]. Note that in these implementations, the viscous term, i.e. the transient term proposed in the present work, is absent. The two PDEs are in fact so similar that in order to solve for the micromorphic variable, instead of the temperature, only two minor modifications need to be done. Given the comparison between the two PDEs (Eqs. 38 and 39), and given the form of the yield function in Eq. 32, the elements that require non-trivial modifications are the source term r and the yield radius: the former has to coincide with the difference between the cumulative plastic strain and the micromorphic variable (amplified by the H parameter), and the latter has to take into account the extra hardening due to the micromorphic variable: whereas the coefficients present in the thermal balance equation can be easily substituted with the parameters characterizing the micromorphic PDE. Implementing the conditions Eqs. 40 and 41 represents the only real, yet minor, effort that is required to make use of the present theory, assuming the existence of a thermal solver and the possibility of applying small modifications.

Influence on the C parameter
The additional parameter C naturally arises from the development of the chosen constitutive material model for the generalized stress a. In order to obtain the final form of the governing Eq. 30, so that the thermal-micromorphic analogy is valid, the presence of the C parameter is required, and it should not vanish in the case of the implementation of the transient problem. However, from the analysis of Eq. 28, it is clear that the parameter C regulates the development of the viscous part of the micromorphic variable, and therefore that a viscous part of the micromorphic variable exists. Being this an additional material parameter, the question on the calibration of such value must be addressed.
The purpose of using the micromorphic analysis, in the present investigation, is to gain indirect control on the distribution of the cumulative plastic strain and its gradient, thus the constraint on the micromorphic variable to closely follow the value of the cumulative plastic strain through the penalty parameter. The present theory also accounts for the development of viscous stresses generated by notnegligible strain rates, and the micromorphic variable follows the value of the cumulative plastic strain, regardless of whether the plastic strain increment is caused by quasistatic or viscous stresses. The adoption of large values of the C parameters (compared to H ) would allow the viscous part of the micromorphic variable to produce additional meaningful generalized stress (see Eq. 28), therefore altering the value that it should have, based only on the difference between micromorphic variable and cumulative plastic strain (effectively producing the same stress as if this difference was larger). Therefore, too large values of C would somehow corrupt and interfere with the equivalence between cumulative plastic strain and micromorphic variable. On the contrary, by neglecting any meaningful contribution of the viscous micromorphic term to exist, we lose the analogy with transient thermal analysis proposed here for the implementation.
Therefore, for the present investigation, the C parameter must exist, so that the thermal-micromorphic analogy holds, but its value should not be too large. The allowed magnitude for this parameter will be tested by checking an analytical solution in the static case, considered in Section 3.1.

Numerical implementation
The micromorphic plasticity model has been implemented in VPS Explicit [10], a finite element software developed by ESI Group solving both dynamics and heat problems. In order to account for the large deformation expected during manufacturing operations, the theory has been developed according to the VPS standard method, that is, using rate-type constitutive equations. This does not alter the theory so far presented, since the micromorphic part remains unchanged. For the same reason, the space gradients that are encountered in this manuscript are meant to be evaluated with respect to the current configuration of the medium, as in an Updated-Lagrangian approach. The additive decomposition is applied to the strain rate tensor, which can be split into elastic and plastic contributions: where D ∼ is the strain rate, and the elastic constitutive model is rewritten by means of a hypoelasticity relation: is the Jaumann stress rate, and it can be re-written as: where W ∼ is the spin tensor. The finite element solution is obtained by establishing the weak form of Eqs. 11 and 12 using the Galerkin method. The dynamic balance Eq. 11 is weighted with the test velocities u whereas the micromorphic balance Eq. 12 is weighted with the test micromorphic variable rates ṗ . Integration over the domain is achieved by the use of the divergence theorem to lower the order of the derivatives. The natural boundary conditions are incorporated as forcing terms, leading to the equations to be discretized by finite-element interpolations. The discretization of the displacement and micromorphic fields over the domain is achieved by using proper-order interpolation functions. The following algebraic equations are derived: where M ∼ is the mass matrix, F ext is the vector of external nodal forces, F int is the vector of internal nodal forces, C ∼ is the viscosity parameter matrix, a r is the vector containing the nodal generalized forces generated by the source terms and a int is the vector of nodal generalized forces induced by Laplacian of the micromorphic variable. In Eq. 46 the similarity with the discretized algebraic equation to solve the heat equation in thermal analysis can be appreciated once again. In fact, VPS Explicit uses the same form of equation to solve the heat equation: where is the nodal temperature vector, C ∼ is the heat capacity matrix, Q is the nodal heat flow depending on the heat flux on the outer surface , Q is the nodal heat flow depending on the internal heat source and Q K is the internal nodal heat flow depending on the heat flux inside the domain .
A central difference explicit scheme associated to the lumped mass matrix is used to solve Eq. 45. Assuming that the problem is initially found at time t 0 and that the objective is to evaluate its status at time t 1 = t 0 + Δt , the following standard steps are taken: A forward Euler scheme associated with the viscosity lumped matrix is implemented to solve Eq. 46: A weak micromorphic-mechanical coupling is implemented in VPS Explicit, that is, the two equations are solved separately so that displacements are considered constants while solving for p and vice-versa. The micromorphic field influences the plastic behavior of the continuum (through condition Eq. 40), and, in return, the cumulative plastic strain (the difference between the cumulative plastic strain and the micromorphic variable) acts as a source term in the micromorphic balance equation (in condition Eq. 41).
Regarding the mechanical behavior, a user material routine implements the mechanical model as previously defined. The values of the micromorphic variables at the Gauss quadrature points are interpolated by mean of the interpolation functions from the nodal values. So the user material routine not only integrates the mechanical behavior but also computes the source term H (p − p ) at the Gauss points. Regarding the micromorphic treatment, a specific function is developed inside the thermal solver in order to recover the source term from the material computations previously evaluated. The main algorithmic steps of the explicit resolution over the time increment Δt may be summarized by the following scheme: The critical time step for the time integration of the equations is taken as the minimum between the critical time step for the mechanical and the micromorphic integration. The critical time step for the mechanical integration follows the standard definition, whereas the critical time step for the integration of the micromorphic variable is the the same used for the thermal variable (following the micromorphic-thermal analogy exploited in Section 2.5), but using the micromorphic parameters: where V is the volume of the element and S max is the largest surface of the element.

Analytical solution
The analytical solution is developed for the rate-independent static case as a reference for validation of the FE scheme at the static limit. It is inspired from similar solution proposed by [26,38,39]. Consider a periodic strip made of a thick rectangular plate of the width W along X 1 direction, the length L along X 2 direction, and the thickness T along X 3 direction (Fig. 1) undergoing simple shear. A macroscopic deformation is applied such that where is the periodic displacement fluctuation. Due to equilibrium conditions, the shear stress component is homogeneous so that the equivalent stress eq is invariant along X 1 , X 2 and X 3 , hence (55) eq (X 1 , X 2 , X 3 ) = eq . The yield condition including the linear softening term and the micromorphic contribution (with C = 0 here) can be written as The PDE governing the micromorphic variable is given by Elimination of the variable p in the previous equation by means of the yield condition Eq. 56 leads to the following form of the PDE to be solved for p : In case of linear softening Eq. 58 takes the form where is the characteristic width of the deformation zone. The expressions for constants and can be found in Appendix 7. The PDE Eq. 59 governing p is only valid in the region X 2 ∈ [ − 2 , 2 ] and the solution is of the form For symmetry reasons, p (X 2 ) = p (−X 2 ) leads to 2 = 0 . At the elastic/plastic interfaces, i.e at X 2 = ± 2 , continuity of micromorphic variable p and of the generalized stress normal to the interface M ⋅ X 2 must hold, hence where we make the approximation that p is sufficiently close to p, i.e. that the penalty coefficient is large enough. Combining Eqs. 61 and 62 with 60 leads to Moreover, the equivalent stress is expressed as where is the elastic shear modulus. From the yield condition, p can be replaced by eq − 0 +H p H p +H in Eq. 64 and integra- tion gives an expression for eq as a function of applied macroscopic shear ̄1 2 and then the uniform shear stress writes where is the shear modulus of the material.

FE solution
The FE simulations are performed with a periodic strip. The associated 2D coordinate system and geometry are shown in Fig. 1. The strip has been meshed with 3D 8-nodes elements onto which plane strain conditions were applied by imposing zero out-of-plane displacement to all the nodes. The nodes at the bottom of the strip ( X 2 = −L∕2 ) were clamped along X 1 and X 2 . The nodes on the top surface ( X 2 = L∕2 ) were clamped along X 2 and a Dirichlet type of boundary condition was applied along X 1 whereas the displacements along X 2 were fixed. Linear shape functions have been used to interpolate the nodal fields, and full integration schemes have been used for the material behavior. Numerically, in order to trigger the strain localization in a periodic strip, a small defect is introduced at the centre (Fig. 1). The defect is one element having an initial yield stress 3% less than the matrix. Isotropic elasticity is considered. The material parameters used for the analytical solution and FE simulations are presented in Table 2. Figure 2a and b show the cumulative plastic strain fields with the classical and the micromorphic models using two different mesh discretizations, one coarse and one fine mesh (using 0.010 mm and 0.003 mm thick elements respectively). Given the fact that in this example the shear band is already known to have a thickness of ≈ 0.08 mm, the chosen mesh size lies well within a safe regime for the gradient effects to be properly captured. The classical plasticity model exhibits pathological mesh dependency and width of the shear band always collapse to one element irrespective of the mesh size. In contrast, the width of the formed shear band with the micromorphic model is finite and independent of the mesh size. This indicates the capabilities of the implemented micromorphic theory in an explicit scheme to solve the shear strain localization problem.
Furthermore, the cumulative plastic strain variation along X 2 obtained from the FE solution is validated against the analytical solution developed for the rate-independent case (cf. Eq. 60), see Fig. 3. The FE simulation is validated for ̄1 2 = 0.01 . Moreover, simulations are performed by changing the simulation time while keeping the same applied total shear strains. Fig. 3b shows that the perfect agreement with an analytical solution is obtained for t = 10 sec. which corresponds to low enough strain rate to make the viscous contribution in Eq. 37 negligible. Larger strain rates are observed to limit the localization since the maximum strain in the band decreases for increasing strain rates. Since the total strain is imposed, this means that a higher elastic strain compensates the lower plastic strain which means that stress values are higher.
In order to retrieve the quasi-static solution, also the viscous parameter C had to be chosen small enough. The reason is to minimize as much as possible any viscous-like component of the generalized stress a in Eq. 28 to retrieve the rate-independent solution.
Metals at high temperatures are known to be strain rate sensitive. This effect is generally taken into account by means of an appropriate viscoplastic flow rule, for instance based on a Norton power law. In the present work, rate-independent plasticity only has been considered but the generalization to viscoplasticity is straightforward in the proposed framework. Note that the proposed model presents an additional strain rate sensitivity, via the viscosity parameter C . This will require appropriate calibration for instance using strain field measurements during localization.

Numerical examples
In this section, the applicability of the implemented scalar micromorphic strain gradient theory is tested for two additional cases: a shearing operation process and a bending test. The aim of this section is to exploit the analogy explained in Section 2.5, whose numerical implementation has been previously presented, to prove that simulations of manufacturing operations using the micromorphic continuum under an explicit integration scheme can be successfully performed.
Industry best practice discourages the employment of complex numerical methods to produce simulations, mainly to guarantee a high degree of reliability of the results and computational efficiency in terms of CPU time. Regarding this reasonable concerns, the results that will be presented here are to be considered as proof of the simplicity of the method, which requires only one additional parameter to be calibrated, that is A (see the discussion on the C parameter in Section 2.6). (a) Equivalent plastic strain distribution obtained for a simulation of 10 seconds for two different mesh discretizations (101 vs. 303 elements); (b) equivalent plastic strain distribution obtained with a fine mesh for different total time As previously explained in the introduction, the relevance of the application of regularization procedures in manufacturing operations is vital, especially in cases in which the thermal power has a major presence. Thermal softening can take place when high rates of plastic strain are produced, and similar softening can be reproduced by assigning a negative slope to the hardening function in Eq. 36. The regularization potential of the proposed method is investigated in the shearing operation section. Moreover, one of the missing features of the classical continuum mechanics is the capability of predicting any size effect. This becomes of major relevance whenever the deformation localizes is small regions or in the case of forming of micro-components [21,46]. The ability of the proposed method to capture the size effect is proven in the bending section.

Shearing operation
The shear band formation is a commonly observed phenomenon in manufacturing operations in case of heavy deformation, for instance, high-speed shaping, forging, machining, and several other processes [4,27]. Numerically, the shear band simulation shows spurious mesh dependency when we consider a classical plasticity approach with strain softening. Dynamics combined with viscosity or/and heat conduction are known to provide regularization but the involved length scales are often too small for efficient FE modeling so that strain gradient or micromorphic plasticity is still useful to introduce physically more realistic length scales [41,45].
Shearing operation is most commonly used in the metal forming industries for sheet metal cutting. In this section, the implemented micromorphic approach is used for the regularization of shear band formation in shearing operation.
The shearing operation is performed on a sheet of 5 mm thickness under plane strain conditions with one element across the width. The geometry is shown in Fig. 4. The sheet has been meshed with 3D 8-nodes elements with linear shape functions and full integration schemes. The lower tool is fixed, while velocity is applied to the upper tool in the downward direction. At the initial deformation stage, a linearly increasing velocity up to 4 mm∕sec. is applied. Once the velocity of 4 mm∕sec. is achieved, it is kept constant in the later stage of the deformation. The contact between the deformable sheet and tools is taken into account using a constant coefficient of friction 0.3. The tools are considered as rigid bodies, while the sheet is assigned with an elastoplastic material behavior using linear strain softening. Isotropic elasticity is considered. The used material parameters in the numerical simulations are presented in Table 3.
At first, simulations are performed with classical plasticity using three different mesh discretizations. The elements size of the different meshes in the region of interest span from 0.2 to 0.04 mm respectively, and, as it will be possible to appreciate in Fig. 6, this element size is extremely small if compared with the shear band thickness, thus ensuring that the gradient effect are correctly captured. The limitation of the classical plasticity model, known as pathological mesh dependency in the strain localization problem can be observed from   and e by the contours of the cumulative plastic strain. The magnitude of the cumulative plastic strain is different for two different mesh discretizations, and it increases with finer mesh. Furthermore, the observed width of the shear band is different for two different mesh discretizations and it always collapses to one element size irrespective of the mesh size. In contrast, the formed width of the shear band using the micromorphic approach is finite and does not depend on the mesh density as seen from Fig. 5b and f. In addition, the magnitude of the cumulative plastic strain reaches asymptotic values while reducing the mesh size. Furthermore, the effect of the diffusivity coefficient A on the shear band widths is investigated. Figure 6 shows the variation of cumulative plastic strain for three different values of the gradient parameters A, 128 N, 320 N, and 800 N. As the value of A increases the intensity of plastic strain gradient within the shear region reduces. As expected from the analytical expression for the length scale in Eq. 66, the width of the shear band increases with an increase in the A value. For the three different values of the A parameter, 128 N, 320 N and 800 N, the observed widths of the shear bands are 2.4 mm, 2.8 mm, and 3.5 mm, respectively. If the characteristic length was evaluated through Eq. 66, the values of 3 mm, 5 mm and 8 mm would be the results. The divergence of these values is due to the fact that the boundary conditions are not the same, thus the deformation state is not pure shear.

Bending
The bending test is used to verify that the implemented micromorphic model is able to capture the size effect for hardening plasticity. It is possible to find in literature many studies that experimentally highlighted the presence of extra hardening in the bending moment, whenever the specimen geometry was reaching sub-millimeters dimension, approaching grain size. In 1994 Fleck and co-workers [11] reported hardening behavior in a copper wire under torsion for wire diameters in the order of 10 − 100 micro-meters, whereas tensile tests performed on the same wires found no evidence of size effect. Stölken and Evans [42] designed a micro-bend test to measure the plastic characteristic length scale associated with the strain gradient, subsequently reporting the results pertaining to thin ( 12.5 m ↦ 50 m ) Nickel foils.
In Fig. 7, the geometry and the boundary conditions of the specimen are reported. The specimen has been discretized using 3D type of elements under plane strain conditions. Linear shape functions are used to interpolate nodal values, and full integration scheme is used for the elements. One element spans the 1 mm width and 10 elements span half the thickness of the beam (mesh size of 0.1 mm) so that there should be enough elements to capture the size effect. The left face of the beam is clamped, whereas a material rotation is enforced on the nodes of the right face through a coupling involving the nodes of the right face and an auxiliary node. The resultant bending moment is probed at the auxiliary node. A total rotation of 45 • is applied.
The size effect can be experimentally encountered whenever the geometry of the specimen reduces down to approximately the grain size of the metal. Virtually, the same phenomenon could be achieved by keeping constant the geometry of the specimen and simultaneously increasing the characteristic length scale. The effectiveness of the formulation in predicting the size effect through the bending test has been verified by employing the latter method. The numerical framework previously presented does not explicitly make use of the grain size, but a characteristic length scale in Eq. 31 was identified, and this will serve the same purpose. The use of larger or smaller characteristic lengths will respectively induce a stiffer or softer global response of the specimen. Three different values of gradient parameter A have been used. The other material parameters used in the simulation of the bending tests are reported in Table 4.  In the attempt of replicating a quasi-static bending test, the chosen value of the C parameters is relatively small, so that any viscous contribution of the micromorphic variable would be negligible. In Fig. 8a the distribution of the cumulative plastic strain for the bending test using the micromorphic model is reported. Besides the edge effect induced by the boundary condition at the right surface, the solution appears to be invariant along the longitudinal direction of the strip.
The classical and micromorphic solutions in terms of normalized bending moment vs. applied rotation are shown in Fig. 8b. The probed bending moment has been normalized with respect to the first moment of area of the beam cross-section, that is b h 2 , where b is the width of the rectangular cross-section, and h is the height of the rectangular cross-section. From Fig. 8b, it can be appreciated that the classical solution is retrieved by using the micromorphic approach with a null penalty term H and null higher-order modulus (A). Three values of the higherorder modulus (respectively three different characteristic lengths scales) are used for the test: 128 N, 320 N, and 800 N. The curves belonging to the micromorphic theory clearly demonstrate the ability of the method to capture the size effect. The extra hardening reported in Fig. 8b follows the same trend as the one relative to the experimental tests reported by Stölken and Evans [42].
In the case of bending, the micromorphic medium does not need to regularize any localization phenomenon; rather, it has to predict an additional hardening, as presented in the manuscript. The characteristic length scale can be identified in this case by l ch = √ A(H p +H ) |H p |H . The obtained characteristic length scales using A = 128 N, 320 N, and 800 N are 0.8 mm, 1.26 mm, and 2.0 mm respectively. The chosen mesh size ensures that many elements are contained within the characteristic length. These characteristic length scales can be normalized by the thickness h of the beam. The obtained l ch ∕h ratios for A = 128 N, 320 N and 800 N are 0.40, 0.63 and 1.0, respectively. Figure 8b shows that for high l ch ∕h ratio, i.e. high A value, stronger response can be predicted. The plasticity material model used for the bending test is characterized by a linear hardening behavior (Table 4). From the analysis of the curves, it can be inferred that the regularization, and subsequently the size effect, is affecting the solution only in the plastic regime, whereas the initial elastic stiffness of the curves is the same regardless of the characteristic length scale used in the model. This is the expected behavior, given the fact that the present micromorphic theory regulates the localization of the plastic field. Thus there should be no difference between the curves in the elastic regime. In hardening plasticity, the plastic strain gradient contribution leads to an increased apparent hardening of the beam in the plastic regime. The peak in terms of adimentionalized bending moment in Fig. 8b is to be attributed to the dynamic nature of the test, and it affects all the curves regardless of the classical or micromorphic nature of the theory used.

Conclusion
In this manuscript, a micromorphic strain gradient plasticity model has been formulated and implemented in a commercial explicit finite element code in order to perform simulations of manufacturing operations in timedependent environments. The reasons to account for the strain gradient while simulating manufacturing operations deal with regularization of strain localization phenomena in softening plasticity, on the one hand, and prediction of size effects in hardening plasticity. The originality of the approach lies in the use of the micromorphic model instead of strict strain gradient plasticity and in the introduction of a viscosity contribution to the micromorphic plastic evolution. The advantage of these two ingredients is that they ease the numerical implementation in a commercial finite element code by mimicking the transient heat equations. Earlier formulations are based on strict strain gradient plasticity without transient term, on the one hand, or on the introduction of micromorphic inertia instead of the proposed viscous term.
The main outcome of the present research lies in the proof that it is possible to implement an explicit micromorphic model in a relatively easy and straightforward manner. This was achieved by slightly modifying the preexisting routines of material integration and thermal field resolution in the VPC/PAMCRASH software developed by ESI. This proof of concept is meant to demonstrate that limited effort is required to implement the micromorphic theory in any other software that allows for minor modification in their procedures.
The implemented theory has been demonstrated to recover the analytical solution for a semi-infinite glide layer under quasi-static loading conditions. The supplementary shearing tests highlighted the need to use of the strain gradient theory in case deformation localizes, and the typical extra hardening in bending has also been modeled.
Most importantly, it has been proven that the size effect can be predicted with this method and that manufacturing operations can be simulated with such theory with a limited increase in computational cost and only one additional material parameter (the characteristic length). The same model can therefore be used to address regularization issues in softening plasticity and "smaller is harder" size effects in microforming. Further work should be dedicated to develop case studies involving real material data and more complex 3D specimen geometries. In particular the consideration of adiabatic shear banding can be included in the approach in a way similar to the work done in [33] whereas full coupling with heat conduction phenomenon would require more intrusive programming in the considered commercial code.