A Review on Strain Gradient Plasticity Approaches in Simulation of Manufacturing Processes

Abstract

Predicting the performances of a manufactured part is extremely important, especially for industries in which there is almost no room for uncertainties, such as aeronautical or automotive. Simulations performed by means of numerical methods such as Finite Element Methods represent a powerful instrument in achieving high level of predictability. However, some particular combinations of manufactured materials and manufacturing processes might lead to unfavorable conditions in which the classical mathematical models used to predict the behavior of the continuum are not anymore able to deliver predictions that are in good agreement with experimental evidence. Since the first evidences of the shortcomings of the classical model were highlighted, many non-classical continuum mechanics theories have been developed, and most of them introduce dependencies at different levels with the Plastic Strain Gradient. This manuscript aims at gathering the milestone contributions among the Strain Gradient Plasticity Theories developed so far, with the object of exploring the way they interface with the requirements posed by the challenges in simulating manufacturing operations. Finally, the most relevant examples of the applications of Strain Gradient Plasticity Theories for manufacturing simulations have been reported from literature.

1. Introduction

The topic of Strain Gradient Theories (SGT) experienced a fervent increase of interest from the scientific community in recent years. Although multiple reasons made the SGT one of the hot topics in the field of Continuum Mechanics (CM), the main causes can be identified in the ever-increasing computational power of numerical calculators (thus enabling the implementation of techniques beforehand impossible to use) and in the demand of more sophisticated CM theories which would better predict the medium behavior under specific material deformation conditions. SGTs could be employed whenever the Classical CM framework is not delivering anymore a proper medium behavior description if compared with experimental observations. The features embodied in the SGTs can be of crucial importance in several situations, and in this manuscript we assess the role that they might play for manufacturing process simulations. However, it must be noted that the adoption of such theories increases the computational cost involved in simulating the observed phenomenon.

The most general feature, shared by all the SGTs, is the introduction of one or more scale lengths in the model used to describe the medium behavior. This is done through the definition of additional deformation measures that are meant to capture specific phenomena which cannot be predicted by the classical CM models. In the vast majority of the cases, the magnitude and the influence of these additional deformation measures are negligible if compared with the standard one, i.e., strain. In specific situations, however, the material behavior demonstrates noticeable susceptibility to these supplementary deformation measures, and proper models are required to capture such behaviors through simulations. The gradient of the strain has been widely and unanimously recognized as one of the main feature which needs to be incorporated in the Classical CM theory to predict such behaviors. Among all the SGTs, those that make use of the gradient of the plastic part of the strain are called Strain Gradient Plasticity Theories (SGPT).

Experimental observations show that many metal manufacturing processes are affected by these phenomena, demonstrating noticeable sensitivity toward non-conventional deformation measures [1,2,3,4,5,6,7,8,9,10,11]. The manufacturing techniques in which the strain gradient plays a bigger role result to be those in which the loads are highly localized, e.g., micro-cutting, milling, micro-bending, thin wire drawing, machining, Friction Stir Welding (FSW), sheet stamping, inverse drawing, adiabatic cutting. These machining processes favor the development of rather peculiar conditions in the processed material, where a combination of separate factors contribute to produce a relatively complex scenario.

The present research aims at addressing the following research question: What is the state of the art in using enhanced continuum mechanics and, more specifically, SGPTs, to simulate manufacturing simulations. This research question can be broken down in the following sub-questions:

• What is an enhanced continuum mechanics model?
• Why such models are required for the simulation of manufacturing operations?
• Which enhanced models of continuum mechanics can be used to simulate manufacturing operations?
• Which contributions can be found in literature that already used enhanced model to simulate manufacturing operations?

In Section 2, this multitude of events simultaneously taking place during these processes are singularly assessed, and then the features that represent them are traduced into requirements of a continuum mechanics theory meant to correctly describe them. Section 3 presents the milestones in the development of different SGPTs and in Section 4 the investigations that employed the SGPTs to simulate machining operations will be thoroughly reviewed. Finally, the conclusions of the present review manuscript will present advantages and disadvantages of using different SGPTs to simulate different machining processes.

2. Challenges in Manufacturing Processes Simulations

Among the different metal manufacturing techniques, the ones that require a proper theory to be described are those that modify the material through strain localization, or severe plastic strain development or material removal, or a combination of those. During machining, for example, material is removed from the work piece by inducing large amounts of shear in localized areas called Shear Bands (SBs) (see Figure 1). Due to the high speed of the process (cutting velocities up to 90 m/s) [12,13,14,15], high strain rates (${10}^5$ s${}^{-1}$) are induced in the SBs [16,17,18]. The high levels of stresses and strains give rise to plastic deformations, therefore inducing softening in the material at this same location. In addition to this, since the SBs usually span an extremely narrow area if compared to the global scale of the work piece, the continuum at this location experiences high stress and high strain gradients as well. To complete the picture, it must be mentioned that plastic deformations usually give rise to heat productions, and considering the case in which the rate of heat production is much larger than the heat flow within the material (especially in metals characterized by low thermal conductivity), high temperature fields are retained, and the SBs are referred to as Adiabatic Shear Bands (ASBs). In Figure 2, the crystallographic analysis during a FSW process highlights the modification of the grain size and orientation induced during manufacturing. Shear bands have been reported to appear also during a high-temperature compression test on Nickel super-alloy [15].

On top of the already complex scenario, at the location where the applied forces are transferred from the tool to the workpiece, the material experiences complex behaviors as well. High stress and shear gradients develop in this zone, leading to the same conditions which are found within the SB. Furthermore, the mechanism with which the loads are transferred from the workpiece to the material must be accurately assessed as well. The extreme conditions in which the contact must be modeled are actually different from the ones for which the classical contact models have been developed.

Overall, the problem of simulating machining seems to be composed of several sub-problems, some of them interconnected with each others, and some independent from one another. This complex problem can be more easily assessed if divided into separate and independent blocks. These might be listed as:

• Strain localization in a length scale whose order of magnitude is the same as grain size, therefore violating the limit of validity of material homogeneity;
• Mesh-size hypersensitivity when introducing material softening in the behavior law during the analysis;
• Proper material characterization at high temperatures and high strain rates developing during manufacturing processes;
• The complex material behavior at the tool-workpiece contact location.

These phenomena could take place separately or simultaneously. Several approaches can be found in literature that address these problems individually, but their occurrence, taking place all at once, represents a challenging issue to be addressed. Complex problems, such as the one investigated here, have also addressed through model reduction methods [19].

2.1. Strain Localization and Mesh-Size Dependency

Strain localization in a narrow area whose order of magnitude is comparable to the material grain size compels the adoption of a non-classical CM framework that would be able to model the size effect. In literature, several authors highlighted the inadequacy of the classical CM framework whenever it was used to reproduce experimental data in which the size effects were noticeable [20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36]. Overall, the inconsistencies are more pronounced whenever the dimensions of the specimens become comparable to the order of magnitude of the specimen grain size, or, equivalently, when the deformations localize in areas whose size is comparable to the specimen grain size.

Fleck and Hutchinson performed several static torsion and tensile tests on copper wires of different diameters (12–170 $\mathsf{\mu }$m) [20]. Their results, reported in Figure 3, highlight that smaller diameter wires are characterized by a much stiffer torsional response (scaled with the wire radius), although tensile tests performed on the same specimen demonstrate a very small, thus negligible, size dependency. Similar tests were performed by Liu et al. who reached the same conclusion [37,38]. Fleck and Hutchinson also reported that the size-effect is present in a Vickers micro-indentation test conducted on single crystal tungsten specimens [20]. The performed tests demonstrated a strong size-dependency, as the material hardness doubles by using an indent whose diagonal is one order of magnitude smaller. Similar size-effects in a micro-indentation test have also been reported by many other authors [23,24,25,26,27].

The size-effect has also been experimentally investigated by Stölken and Evans [21]. They presented the results of micro-bending tests conducted on 12.5 $\mathsf{\mu }$m, 25 $\mathsf{\mu }$m and 50 $\mathsf{\mu }$m thin Nickel foils (Figure 4). The results of their tests highlighted that foils with a smaller thickness were behaving stronger. The applied bending moment (normalized with respect to the foil thickness) of the 50 $\mathsf{\mu }$m thick foil is twice the one recorded with a 12.5 $\mathsf{\mu }$m thick foil. Similar tests leading to similar conclusions were conducted by Ehrler et al. [39].

In the same referenced papers, so as in a broad branch of literature [28,29,30,31,32,33,34,35,36], the authors point out that the classical CM is not able to capture these localization phenomena due to the absence of a length scale in the model that would counteract the localization of the fields. Micro-torsion tests also highlighted the shortcoming of the classical CM framework [40,41]. The SGPT, whose mathematical description includes one or more length scales, has been proposed as a valid candidate to capture phenomena of different natures.

The size-effect can be predicted thanks to the definition of an additional deformation measure, i.e., the strain gradient. The dimensional analysis of this quantity reveals that it is not dimensionless and, unlike strain, it should be adimensionalized by means of a specific length, and the sensitivity of non-classical CM theories to size-effects comes directly from the definition of this non-dimensionless deformation measure. Many researchers investigated over the physical nature of the length scales of the SGPTs. Depending on the degrees of complexity of the observed phenomena, more than one characteristic lengths can be identified. For example, Fleck and Hutchinson identified up to five distinct characteristic lengths, each associated to one of the five invariants of the strain gradient tensor [20]. These characteristic lengths however can be reduced to three in case of incompressible isotropic medium. In a successive study, Fleck and Hutchinson identified a single characteristic length in an equivalent plastic strain gradient theory [42]. Liu and Dunstan, based on the SGPT of Fleck and Hutchinson, gave a physical interpretation to the characteristic length by making a connection to physical quantities via critical thickness theory [33]. Duan et al., based on the same SGPT, associated the additional characteristic length to geometrically necessary dislocations through three different dislocation models [35]. Similarly, Dahlberg and Boåsen assumed an equivalence between the microstructural length scale and dislocation density and provided an evolution law of the characteristic length [43]. Zhang and Aifantis gave a comprehensive review of the interpretations of the characteristic length associated with SGPTs [44].

Several other researchers focused on the calibration procedure and the quantification of this length for metals. Yuan and Chen proposed to identify the characteristic length from micro- and nano- indentation tests [45]. Stölken and Evans developed a micro-bending test to measure the characteristic length [21]. Abu Al-Rub and Voyiadjis also proposed to adopt micro- and nano- indentation tests to calibrate the characteristic length and its evolution law [34,46].

Thereby, it can be concluded that the problem of properly capturing the size-effect rising during manufacturing simulation due to strain localization can be solved by adopting a medium description that includes the gradient of the strain as a deformation measure. In addition, adding the strain gradient as a deformation measure would solve not only the problem related to size-effect, but also the one related to mesh hypersensitivity.

Mesh hypersensitivity is an issue which has its root in the loss of ellipticity of the partial differential equation governing the medium equilibrium [47,48,49]. In case the non-linear plastic behavior of the material experiences a decreasing flow stress, i.e., softening, the nature of the set of partial differential equations governing the equilibrium changes, and a strong mesh sensitivity is experienced, in particular, the solution does not appear to converge toward an asymptotic value when the mesh size is decreased. The process of recovering the property of mesh independence in literature is referred to as regularization procedure.

Jirasek and Rolshoven provided an extensive comparison of the regularization properties of many SGPT by analyzing the response of a mono-dimensional bar under tension, whose plastic behavior is characterized purely by softening [50]. They explored the regularization mechanics of the Toupin-Mindlin elastic SGT [51,52], the plastic version of Toupin-Mindlin SGT developed by Chambon et al. [53], the Fleck and Hutchinson plastic SGT [20] and the Mechanism-based SGT [54]. Nguyen et al. coupled a non-local plasticity model with damage to successfully capture the softening behavior experienced by the material during ductile failure in the post-peak regime [55]. Lele and Anand demonstrated that SGPTs are able to provide mesh independent results in case of viscoplastic material as well [56].

Several other researchers overcame the mesh hypersensitivity issue through other solutions besides the SGT. Mediavilla et al. used a damage-enriched material model, in which the gradient of the damage field would enter in the material model, thereby achieving the same regularization effect produced by the SGTs [57]. Many other researchers included the gradient of the damage variable in the model to regularize the solution [55,58,59,60]. Higher Order continuum descriptions (see Section 3) are also well known equivalent solution to avoid the pathological mesh dependence [61,62,63,64]. In Figure 5, proof of the regularization property of a Higher Order theory (Cosserat medium) is reported. A special combination of specimen geometry and applied boundary conditions induce the appearance of SBs inside the material [65], and it has been demonstrated that the localization of the equivalent plastic strain inside the SB is mesh-independent in case a Generalized CM (Cosserat) theory is used.

Therefore, the issues of mesh hypersensitivity due to localization during machining simulations can be simultaneously addressed through the adoption of non-local theories, to which the SGPT belong.

2.2. Material Characterization at High Temperatures and High Strain Rates

The high strain rate fields produced during machining are inevitably intertwined with thermal-related consequences. Materials experiencing high strain rate are likely to plasticize and subsequently to generate heat. If the rate of heat production is higher than the rate of heat conduction, high temperatures will be retained in the plasticized zone. The consequences of high temperatures on material behavior can be easily modeled and they might be independently analyzed and simply assessed by performing material tests at high temperatures; on the contrary, the reliability of the model relating high strain rates with mechanical behavior, due to unavoidable raises in temperature, is dependent on the fidelity of the research done on the thermal behaviors.

In general, it must be recognized that reproducing the same level of shear deformation ($\approx 2.0$) and strain rate ($\approx {10}^5{s}^{-1}$) that takes place during manufacturing processing under controlled lab conditions is not an easy task. The thermal aspect involved in the development of high strain rates was reported by Marchand and Duffy, who collected experimental data of both temperatures and shear strain fields during a dynamic torsion test [66]. They tested a thin hollow tube of HY-100 steel in a torsional Kolsky (split-Hopkinson) bar properly modified to provide torsional loading at high strain rate (1600 s${}^{-1}$). The specimen was 0.38 mm thick and 2.5 mm long, with an inner diameter of 9.5 mm. The authors adopted such high strain rates so as to induce the formation of adiabatic shear bands and then measure the temperatures that developed in the ASB. This condition is equivalent to the one that can occur during machining processes. From the collected data, reported in Figure 6, it can be inferred that the thermal contribution to the energy balance cannot be neglected when SBs develop during machining.

The influence of the thermal field and strain rates on the mechanical behavior was also recently investigated by Rahmaan et al. [67] (Figure 7). They performed dynamic shear tests on 2 mm thick AA7075-T6 sheets at different strain rates (from 0.01 to 10${}^3$ s${}^{-1}$) inducing increments in temperatures of up to 60 °C in the specimen due to plasticity. The authors distinguished between a strain hardening dominant region and thermal softening dominant region. From their results, it is evident that plastic strain at high strain rates (500 s${}^{-1}$) induce a temperature raise causing material softening. Additionally, the temperatures profiles obtained using different strain rates were directly compared, and the analysis showed that the influence of the strain rate on the developed temperature is considerably lower when the strain rates are reduced.

The researches reported in this section proves the fact that in the sought of a complete theory to model machining, the thermal and dynamic contributions cannot absolutely be neglected. Therefore, besides the usual governing equations relative to force equilibrium, a thermodynamical approach must be adopted so that energetic equilibrium equations would be included in the CM model. Therefore, the simulations of manufacturing processes involving very high deformation rates require an approach that would encompass thermodynamical considerations. This can be achieved by properly characterizing the media behavior from a deeper thermodynamical standpoint (from which the material constitutive behavior naturally derives), while solving the energy equation (heat equation) in addition to the standard set of equations governing the force equilibrium. This approach requires the temperature to be included as an additional degree of freedom of the continuum whose field is the solution of the energy equation. However, if the heat dissipation rate is very small compared to the head production rate, the boundary value problem can be considered as adiabatic, thus demoting the temperature from a degree of freedom of the system to a state variable [9,13,16,68].

Another crucial consideration is that thermal-induced effects could be isolated and separately analyzed, but the opposite cannot be done with high strain rates, which are unavoidably related with heat productions. High temperatures soften the material, whereas higher strain rate have an opposite hardening effect. However, if high levels of strains are rapidly obtained, the material plasticize and temperature developments are expected to occur, thus interconnecting the strain rate effect with the thermal effect. So the ideal approach would be to separate these two aspects: first developing a model of the mechanical behavior that delivers the correct effect of thermal variation on the medium response, and then to successively incorporate the effect of the strain rate on mechanical behavior in terms of hardening.

2.3. Tool-Workpiece Contact

During machining processes, the contacts conditions between the tool and the workpiece lie beyond the assumptions made for the contact models usually used in mechanics. Especially the high temperatures and the high pressures, in fact, induce the material to behave more liquid-like than solid-like (see Figure 8). The simple Coulomb law is not a good model anymore [69].

Zorev’s model is a well-known refined friction model that is widely used for machining simulation. It recognizes two different zones at the tool-workpiece contact zone, each characterized by two different frictional behavior [70]. He individuates a sliding zone, where the material response is purely elastic, and a sticking zone, that is instead characterized by plastic flow.

This special condition does not only require a modification of the classical CM governing equations, but it also impels the adoption of a proper numerical contact model able to cope with a change in the friction conditions from a sliding form to a sticking form (transversely promoted by high temperatures and stresses also). Based on this model, further modifications have been reported in literature attempting at delivering a precise model characterizing the metal behavior in the sliding and sticking regime [71,72,73,74,75].

A throughout review of experimental data on the topic of tool-workpiece interaction was collected by Astakov [69]. In this book chapter he reported many experiments which attempt to extrapolate many features of the phenomena from a tribological point of view such as the friction coefficient (assuming Coulomb law), contact stress distribution, thermal distribution and more.

Towards the development of a more appropriate contact model, the FSW benchmark, as suggested by Cahuc et al. [7], might result in a useful mean which could be used to test the validity of the model itself. However, being FSW characterized by some on the main above-cited typical problems affecting machining, other independent features must be assesses before being able to test a more complex model.

The enhancement of the governing equations of the continuum mechanics from a classical to a generalized framework will likely bring the introduction of more complex variables associated with status of the continuum, such as strain gradients. Therefore, further considerations must be made on how these additional unknowns will fit into the proposed contact model, or even more important, how the contact conditions affect material deformations in enhanced continua.

3. Historical Excursus of SGPTs in Literature

The SGPTs belong to a branch of mechanics that aims to describe the continuum behavior through a set of equations which differ from the one belonging to the Classical CM [76,77]. This class of theories is commonly referred to as Generalized Continuum Mechanics GCM. This definition is attributed to a quite large group of theories, and they can be further grouped into other two categories: Higher Order Theories and Higher Grade Theories [64,78]. The former group aims to enhance the continuum with additional degrees of freedom, and the latter adopts the same quantities characterizing the classical CM framework, in addition to their higher spacial gradients. Although the distinction here is quite net, some proposed theories present features that are common to both classes. Most SGPTs belong to the Higher Order Theories, because they include partial or full contributions of the second gradient of the displacements.

Another major distinction that can be drawn among the GCM theories consists in using the additional deformation measure as internal variable or as additional degree of freedom of the media, and this choice has significant repercussions also at the level of numerical implementation. In the first case, the only degrees of freedom of the continuum are the displacements, and the additional deformation measure is evaluated purely from them. In the second case, the additional deformation measure is used as a degree of freedom (or it derives from the additional degrees of freedom) and its distribution resolves as the solution of the set of partial differential equation, so as the displacement field. The Strain Gradient Theories that follow are here explored in their small deformation framework; tensorial notation is reported in the Appendix A.

3.1. Aifantis’ Theory

The pioneering work of Aifantis belongs to a class of SGTs known as Mechanics-Based Strain Gradient (MSGT). The original idea of Aifantis was to enhance the classical Equation for the flow stress with a term related to the spacial gradient of the plastic strain in order to represent the effect of dislocation-state evolution on a macroscopic scale [79,80]. The yield stress as defined by Aifantis reads [36]:

$${\sigma }^Y={\sigma }_0+k{\epsilon }_{eq}^p-c_1|\nabla {\epsilon }_{eq}^p|-c_2{\nabla }^2{\epsilon }_{eq}^p;$$

(1)

where ${\sigma }^Y$ is the flow stress, ${\sigma }_0$ is the initial flow stress, $k,c_1$ and $c_2$ are parameters, ${\epsilon }_{eq}^p$ is the von Mises cumulative plastic strain and ∇ is the nabla operator, returning the gradient/Laplacian of the field. By using this flow stress evolution law, the strain field in the medium has to be found through solution of partial differential equation. The two characteristic lengths are embedded in the definition of the constants $c_1$ and $c_2$. The work of Aifantis was one of the first researches linking micro-defects, such as dislocations, to macro-deformation measures, such as strain. However, the definition of the internal power of the theory was not explicitly stated, and it was not declared whether the inelastic deformations were dissipative or not (the thermodynamic of the model was not originally defined). Successively, Gudmundson [81] incorporated the Aifantis model into the broader framework that he proposed, identifying the internal power density definition as:

$$p^{\left(i\right)}=\underset{\sim }{\mathsf{\sigma }}:{\dot{\underset{\sim }{\mathsf{\epsilon }}}}^e+q{\dot{\mathsf{\epsilon }}}_{eq}^p+\underline{\mathbf{m}}·\left(\underline{\mathbf{\nabla }}{\mathsf{\epsilon }}_{eq}^p\right);$$

(2)

where the stress has been partitioned into the contribution conjugate to the elastic deformation ($\underset{\sim }{\mathsf{\epsilon }}$), $\underset{\sim }{\mathsf{\sigma }}$, and the part conjugate to the equivalent inelastic strain (${\mathsf{\epsilon }}_{eq}^p$), that is a scalar q, and, in addition, the gradient of the equivalent plastic strain is introduced into the internal power definition, where its power conjugate, the moment stress $\underline{\mathbf{m}}$, naturally arises. Gudmundson interpreted this stress as derived from a free energy potential in a quadratic form, as originally suggested by Steinmann [82]:

$$\psi =L^{2}E\underline{\mathbf{\nabla }}{\mathsf{\epsilon }}_{eq}^p·\underline{\mathbf{\nabla }}{\mathsf{\epsilon }}_{eq}^p/2;$$

(3)

from which, the constitutive model for the moment stress derives as:

$$\underline{\mathbf{m}}=L^{2}E\underline{\mathbf{\nabla }}{\mathsf{\epsilon }}_{eq}^p;$$

(4)

Later on, the thermodynamical considerations related to the Aifantis model have been investigated by Gurtin and Anand [83], and they outlined that the Aifantis model cannot be thermodynamically compatible in absence of a defect energy.

3.2. Gradient of the Local Spin Vector—Fleck and Hutchinson 1993

Fleck and Hutchinson developed a phenomenological SGT based on the couple-stress theory [51,52,84], with the idea of relating the strain gradient to the dislocations developed at microstructural level [85]. Besides the standard definition of strain, they considered also the effects of the skew-symmetric part of the displacement gradient, i.e., the local spin tensor $\underset{\sim }{\mathsf{\theta }}$, defined as:

$$\underset{\sim }{\mathsf{\theta }}=\mathrm{skew}\left(\underline{\mathbf{u}}\otimes \underline{\mathbf{\nabla }}\right);$$

(5)

being this a skew-symmetric tensor, it can rewritten in a pseudo-vectorial form as:

$$\underline{\stackrel{\times }{\mathsf{\theta }}}=-{\displaystyle \frac{1}{2}}\underset{\simeq }{\mathsf{ϵ}}:\underset{\sim }{\mathsf{\theta }};$$

(6)

where $\underset{\simeq }{\mathsf{ϵ}}$ is the Levi-Civita permutation tensor. The gradient of the local spin vector, which could be though as a material curvature gradient, is:

$$\underset{\sim }{\mathsf{\chi }}=\underline{\stackrel{\times }{\mathsf{\theta }}}\otimes \underline{\mathbf{\nabla }};$$

(7)

and an additional decomposition in plastic and elastic part of curvature is proposed:

$$\underset{\sim }{\mathsf{\chi }}={\underset{\sim }{\mathsf{\chi }}}^e+{\underset{\sim }{\mathsf{\chi }}}^p;$$

(8)

Based on the definition of the given deformation measures, the internal power would assume the following shape:

$$p^{\left(i\right)}=\underset{\sim }{\mathsf{\sigma }}:\dot{\underset{\sim }{\mathsf{\epsilon }}}+\underset{\sim }{\mathbf{m}}:\dot{\underset{\sim }{\mathsf{\chi }}};$$

(9)

where $\underset{\sim }{\mathsf{\sigma }}$ is the Cauchy stress, energetically conjugate to the strain, thereby symmetric tensor, and $\underset{\sim }{\mathbf{m}}$ is the not-symmetric couple stress tensor as also introduced by Toupin and Mindlin [51,52,84]. The definition of the strain energy density function for an isotropic material was given in a quadratic form as:

$$u^{\left(i\right)}={\displaystyle \frac{\lambda }{2}}{\epsilon }_{ii}^e{\epsilon }_{jj}^e+\mu {\epsilon }_{ij}^e{\epsilon }_{ij}^e+\mu l_{el}{\chi }_{ij}^e{\chi }_{ij}^e;$$

(10)

where $\lambda$ and $\mu$ are the Lamé parameters and $l_{el}$ is an elastic characteristic length used to a-dimensionalize the curvature tensor. The plastic behavior is governed by an equivalent plastic deformation measure that was defined as:

$${\Xi }^{p^2}={\epsilon }_{eq}^{p^2}+l_{pl}^2{\chi }_{eq}^{p^2};$$

(11)

where ${\epsilon }_{eq}^p$ is the von Mises strain invariant $\sqrt{2/3{\epsilon }_{ij}^p{\epsilon }_{ij}^p}$ and similar definition corresponds to ${\chi }_{eq}^p$. The length $l_{pl}$ is used to a-dimensionalize the plastic contribution of the material curvature tensor (defined as gradient of the plastic part of the local spin tensor), and it was though as the free-slip distance between statistically stored dislocations on an average sense. If the condition $l_{el}\ll l_{pl}$ is applied, then the material curvature is assumed to be mostly plastic, therefore providing a sensitive contribution in characterizing plastic phenomena, while being negligible in elastic regimes.

3.3. Second Gradient of Displacement—Fleck and Hutchinson 1997

In 1997 Fleck and Hutchinson enriched the Mindlin-Toupin theory [51,52,84] with a J2 plasticity contribution [20]. Instead of considering only the gradient of the material spin tensor as done in [85], they used the full second gradient of the displacement vector $\underset{\simeq }{\mathsf{\eta }}=\underline{\mathbf{u}}\otimes \underline{\mathbf{\nabla }}\otimes \underline{\mathbf{\nabla }}$:

$$p^{\left(i\right)}=\underset{\sim }{\mathsf{\sigma }}:\dot{\underset{\sim }{\mathsf{\epsilon }}}+\underset{\simeq }{\mathsf{\tau }}⠇\dot{\underset{\simeq }{\mathsf{\eta }}};$$

(12)

where $\underset{\sim }{\mathsf{\sigma }}$ is the Cauchy stress, $\underset{\sim }{\mathsf{\epsilon }}$ is the strain, $\underset{\simeq }{\mathsf{\tau }}$ is commonly referred to as higher order stress and ⠇ is the triple contraction. By assuming an additive decomposition of the second gradient of displacement:

$$\underset{\simeq }{\mathsf{\eta }}={\underset{\simeq }{\mathsf{\eta }}}^e+{\underset{\simeq }{\mathsf{\eta }}}^p$$

(13)

the definition of the strain energy density function is given once again in a quadratic form:

$$u^{\left(i\right)}={\displaystyle \frac{\lambda }{2}}{\epsilon }_{ii}^e{\epsilon }_{jj}^e+\mu {\epsilon }_{ij}^e{\epsilon }_{ij}^e+a_1{\eta }_{ijj}^e{\eta }_{ikk}^e+a_2{\eta }_{iik}^e{\eta }_{kjj}^e+a_3{\eta }_{iik}^e{\eta }_{jjk}^e+a_4{\eta }_{ijk}^e{\eta }_{ijk}^e+a_5{\eta }_{ijk}^e{\eta }_{kji}^e;$$

(14)

where $a_i$ are material parameters to be calibrated and they include characteristic lengths. For incompressible medium, some of the terms in the internal energy definition would drop, and the internal energy would be function of the terms associated to $\mu ,a_3,a_4$ and $a_5$ only [86]. An equivalent deformation measure is then defined as function of the remaining terms, and it will characterize the standard plastic behavior of the medium:

$${\Xi }^{p^2}={\epsilon }_{eq}^{p^2}+\sum _{i=1}^3{l}_{{pl}_i}^2{\chi }_{{eq}_i}^{p^2};$$

(15)

where $l_{{pl}_i}$ are three additional plastic characteristic lengths in the model, which will govern the three different deformation modes each of one embodied by the three remaining contributions of the second displacement gradient in the equivalent deformation measure definition, that are, ${\chi }_{{eq}_i}^{p^2}$.

3.4. Gradient of the Cumulative Plastic Strain—Fleck and Hutchinson 2001

Successively, in 2001, Fleck and Hutchinson proposed even another version of the theory, in which the additional variable was not the whole second gradient of displacement (that also includes the gradient of the local rotation), but they focused on the gradient of the plastic strain only [42]:

$${\underset{\sim }{\mathsf{\epsilon }}}^p\otimes \underline{\mathbf{\nabla }};$$

(16)

A direct comparison can be made between the SGPTs of Fleck1997 and Fleck2001 if we compare the definition of the additional deformation measures they enriched their theory with:

$${\left[\mathrm{sym}\left(\underline{\mathbf{u}}\otimes \underline{\mathbf{\nabla }}\right)\right]}^p\otimes \underline{\mathbf{\nabla }}\ne {\left(\underline{\mathbf{u}}\otimes \underline{\mathbf{\nabla }}\otimes \underline{\mathbf{\nabla }}\right)}^p;$$

(17)

This comparison highlights that the gradient of the plastic strain does not necessarily coincides with the plastic part of the second gradient of displacement. In this case, the strain energy potential would not depend on the gradient of the plastic strain, effectively achieving the same description as in Fleck1993 if the characteristic length was set to be much smaller than the plastic one:

$$u^{\left(i\right)}={\displaystyle \frac{\lambda }{2}}{\epsilon }_{ii}^e{\epsilon }_{jj}^e+\mu {\epsilon }_{ij}^e{\epsilon }_{ij}^e;$$

(18)

while the total power would still feels the contribution of the gradient:

$$p^{\left(i\right)}=\underset{\sim }{\mathsf{\sigma }}:\dot{\underset{\sim }{\mathsf{\epsilon }}}+\underset{\simeq }{\mathsf{\psi }}⠇\left({\dot{\underset{\sim }{\mathsf{\epsilon }}}}^p\otimes \underline{\mathbf{\nabla }}\right);$$

(19)

where $\underset{\simeq }{\mathsf{\psi }}$ is an higher order stress that activates only in case gradients of plastic deformation are presents. The equivalent plastic strain in the present theory measures:

$${\Xi }^{p^2}={\epsilon }_{eq}^{p^2}+l_{pl}^2{|\underline{\mathbf{\nabla }}{\epsilon }_{eq}^p|}^2;$$

(20)

However, two subsequent contributions from Gurtin and Anand [83] and Gudmundson [81], when analyzing the thermodynamical compatibility of this model, found out that it violates the requirement on the positiveness of the dissipation. Following these considerations, Fleck and Willis first [87], and Hutchinson later [88], provided a reformulation of the theory with the constitutive model of incremental nature, so that it fits in a thermodynamically-compatible framework.

3.5. Irrotational Plastic Flow and Burgers Tensor—Gurtin and Anand 2005

In 2005 Gurtin and Anand proposed a SGPT characterized by an irrotational plastic flow and by Burgers tensor as an energetically deformation measure [32,89]. Being the plastic strain irrotational, it means that, contrary to Fleck1993, the condition ${\underset{\sim }{\mathsf{\theta }}}^p=0$ is valid, therefore this term has no energetic contribution:

$$p^{\left(i\right)}=\underset{\sim }{\mathsf{\sigma }}:{\dot{\underset{\sim }{\mathsf{\epsilon }}}}^e+\underset{\sim }{\mathsf{\tau }}:{\dot{\underset{\sim }{\mathsf{\theta }}}}^e+{\mathsf{\sigma }}^p:{\dot{\underset{\sim }{\mathsf{\epsilon }}}}^p+\underset{\simeq }{\mathsf{\psi }}⠇\left({\dot{\underset{\sim }{\mathsf{\epsilon }}}}^p\otimes \underline{\mathbf{\nabla }}\right);$$

(21)

The definition of the density of free energy of the medium was explicitly made dependent on elastic strain, ${\underset{\sim }{\mathsf{\epsilon }}}^e$, and Burgers tensor, $\underset{\sim }{\mathbf{G}}$, that reads:

$$\underset{\sim }{\mathbf{G}}=\mathrm{curl}{\underset{\sim }{\mathsf{\epsilon }}}^p;$$

(22)

therefore, this means that the part of $\underset{\simeq }{\mathsf{\tau }}$ associated with $\underset{\sim }{\mathbf{G}}$ is recoverable, and it plays the role of a back-stress. They identify a dissipative characteristic length, $l_{pl}$, that comes into plays when an equivalent plastic deformation is defined as:

$${\Xi }^{p^2}=|{\underset{\sim }{\mathsf{\epsilon }}}^p{|}^2+l_{pl}^2{|\left({\underset{\sim }{\mathsf{\epsilon }}}^p\otimes \underline{\mathbf{\nabla }}\right)|}^2;$$

(23)

together with an energetic characteristic length, $l_{el}$, that associates the deformation measures in the definition of the internal recoverable energy in a quadratic form:

$$u^{\left(i\right)}={\displaystyle \frac{\lambda }{2}}{\epsilon }_{ii}^e{\epsilon }_{jj}^e+\mu {\epsilon }_{ij}^e{\epsilon }_{ij}^e+{\displaystyle \frac{1}{2}}\mu l_{el}^2{G}_{ij}{G}_{ij};$$

(24)

3.6. The Common Framework—Gudmundson 2004

Gudmundson provided a general framework with the potential to incorporate many SGPTs [81]. The internal power definition is his model depends of the elastic strain, full plastic strain tensor and the gradient of the full plastic strain:

$$p^{\left(i\right)}=\underset{\sim }{\mathsf{\sigma }}:\dot{\underset{\sim }{\mathsf{\epsilon }}}+\left(\underset{\sim }{\mathbf{q}}-{\underset{\sim }{\mathsf{\sigma }}}^{\prime }\right):{\dot{\underset{\sim }{\mathsf{\epsilon }}}}^p+\underset{\simeq }{\mathsf{\psi }}⠇\left({\dot{\underset{\sim }{\mathsf{\epsilon }}}}^p\otimes \underline{\mathbf{\nabla }}\right);$$

(25)

where $\underset{\sim }{\mathbf{q}}$ is a deviatoric microstress, power conjugate to the plastic strain, ${\underset{\sim }{\mathsf{\sigma }}}^{\prime }$ is the deviatoric part of the Cauchy stress and $\underset{\simeq }{\mathsf{\psi }}$ is an higher order stress, power conjugate to the gradient of the plastic strain. From the application of the divergence theorem, it is possible to write the equilibrium governing equations of the medium:

$$\begin{array}{c}\underset{\sim }{\mathsf{\sigma }}·\underline{\mathbf{\nabla }}=\underline{0};\end{array}$$

(26)

$$\begin{array}{c}\underset{\simeq }{\mathsf{\psi }}·\underline{\mathbf{\nabla }}+{\underset{\sim }{\mathsf{\sigma }}}^{\prime }-\underline{\mathbf{q}}=\underline{0};\end{array}$$

(27)

Additionally, Gudmundson defined the constitutive equations through thermodynamical considerations, by explicitly evaluating the medium dissipation and by applying the Clausius–Duhem inequality. The innovation of the theory lies in the identification of a stress conjugate to the plastic strain and plastic strain gradient from the thermodynamical point of view, such that the dissipation can be easily defined as:

$$\left(\underset{\sim }{\mathsf{\sigma }}-\frac{\partial \varphi }{\partial {\underset{\sim }{\mathsf{\epsilon }}}^e}\right):{\dot{\underset{\sim }{\mathsf{\epsilon }}}}^e+\left(\underset{\sim }{\mathbf{q}}-\frac{\partial \varphi }{\partial {\underset{\sim }{\mathsf{\epsilon }}}^p}\right):{\dot{\underset{\sim }{\mathsf{\epsilon }}}}^p+\left[\underset{\simeq }{\mathsf{\psi }}-\frac{\partial \varphi }{\partial \left({\underset{\sim }{\mathsf{\epsilon }}}^p\otimes \overline{\mathbf{\nabla }}\right)}\right]:\left({\dot{\underset{\sim }{\mathsf{\epsilon }}}}^p\otimes \underline{\mathbf{\nabla }}\right);$$

(28)

where $\varphi$ is the Helmholtz free energy potential. By imposing the second thermodynamic principle (positive dissipation), Equation (28) defines the constitutive model of the medium (elastic part) and the dissipation:

$$\begin{array}{c}\underset{\sim }{\mathsf{\sigma }}=\frac{\partial \varphi }{\partial {\mathsf{\epsilon }}^e};\end{array}$$

(29)

$$\begin{array}{c}\pi =\left(\underset{\sim }{\mathbf{q}}-\frac{\partial \varphi }{\partial {\underset{\sim }{\mathsf{\epsilon }}}^p}\right):{\dot{\underset{\sim }{\mathsf{\epsilon }}}}^p+\left[\underset{\simeq }{\mathsf{\psi }}-\frac{\partial \varphi }{\partial \left({\underset{\sim }{\mathsf{\epsilon }}}^p\otimes \overline{\mathbf{\nabla }}\right)}\right]:\left({\dot{\underset{\sim }{\mathsf{\epsilon }}}}^p\otimes \underline{\mathbf{\nabla }}\right);\end{array}$$

(30)

where $\pi$ is the medium density dissipation. He showed that, with this mathematical description, it is possible to retrieve the Aifantis1987 SGPT, the Fleck1997 SGPT and Fleck2001 SGPT.

3.7. Dislocations-Enriched SGT

Gao and coworkers proposed in 1999 a Mechanism-Based Strain Gradient Theory (MSGT) [90,91]. This SGT aims at correlating the strain gradient with Statistically Stored Dislocations (SSD) and Geometrically Necessary Dislocations (GND), and it is built upon the dislocations-based Taylor’s flow stress [92], whose shear stress law description reads:

$$\sigma =\alpha \mu b\sqrt{{\rho }_{SSD}+{\rho }_{GND}};$$

(31)

where $\sigma$ is the flow stress, $\alpha$ is a constant usually assumed to be 0.5 [24], $\mu$ is the shear modulus and b is the magnitude of the Burgers vector. In a single crystal, the GNDs are directly related with the plastic strain gradient [21,93,94]:

$${\rho }_{GND}={\displaystyle \frac{\overline{r}{\eta }_{eq}}{b}};$$

(32)

where $\overline{r}$ is the Nye factor and ${\eta }_{eq}$ is an equivalent measure of the strain gradient, that they evaluated as:

$${\eta }_{eq}=\sqrt{{\displaystyle \frac{1}{4}}\underset{\simeq }{\mathsf{\eta }}⠇\underset{\simeq }{\mathsf{\eta }}};$$

(33)

The flow stress, assuming von Mises plasticity, can be written as [24]:

$$\mathsf{\sigma }={\mathsf{\sigma }}_Y\sqrt{f^2\left(\epsilon \right)+l{\eta }_{eq}^p};$$

(34)

where ${\mathsf{\sigma }}_Y$ is the yield stress measured without strain gradient effect (${\rho }_{GND}=0$), $f^2\left(\epsilon \right)$ is a function of the strain, and l is a characteristic length that reads:

$$l=3{\alpha }^2{\left({\displaystyle \frac{\mu }{{\mathsf{\sigma }}_Y}}\right)}^{2}b;$$

(35)

The introduction of the strain gradient in the MSGT, entails the introduction of an additional term in the internal power density definition:

$$p^{\left(i\right)}=\underset{\sim }{\mathsf{\sigma }}:\dot{\underset{\sim }{\mathsf{\epsilon }}}+\underset{\simeq }{\mathsf{\psi }}⠇\dot{\underset{\simeq }{\mathsf{\eta }}};$$

(36)

where $\underset{\sim }{\mathsf{\sigma }}$ and $\underset{\sim }{\mathsf{\epsilon }}$ are the usual Cauchy stress and usual strain, and $\underset{\simeq }{\mathsf{\psi }}$ and $\underset{\simeq }{\mathsf{\eta }}$ are higher order stress tensor and strain gradient tensor. One of the conclusions outlined by Gao et al. was that due to the fact that:

$${\displaystyle \frac{\partial {\sigma }_{ij}}{\partial {\eta }_{klm}}}\ne {\displaystyle \frac{\partial {\psi }_{klm}}{\partial {\epsilon }_{ij}}}$$

(37)

a potential for the strain energy density cannot exist. In the MSGT, the size effect is induced by including the influence of the strain gradient in the flow stress model.

In 2004 Huang proposed a simplified model, based on the MSGT, called the Conventional Mechanism-Based Strain Gradient Theory (CMSGT). This theory is called a low-order theory, whereas the MSGT is a high order theory, and the difference lies in the fact that the CMSGT makes use of the plastic strain gradient only at the constitutive level, leaving the contribution of the higher order stress out of the theory. This notably simplifies the implementation of the theory, although the well-posedness of the associated boundary value problem is not ensured [95].

Another major contribution in the development of dislocations-based strain gradient theories can be identified in the work of Menzel and Steinmann [96]. They proposed a phenomenological description of the hardening process through introduction of second and fourth order gradient of the plastic deformation for single crystals and polycrystals materials respectively. They exploited the second thermodynamic principle (Clausius–Duhem inequality), in which the definition of the dissipation was given in an additive decomposition. One of the terms composing the dissipation was dependent upon the dislocation density, which was measured as incompatibility of the plastic part of the gradient of the displacement field:

$${\underset{\sim }{\mathsf{\alpha }}}^p=\mathrm{curl}{\left(\underline{\mathbf{u}}\otimes \underline{\mathbf{\nabla }}\right)}^p;$$

(38)

As an example, they proposed, for single crystals, a quadratic Helmholtz potential function of ${\underset{\sim }{\mathsf{\alpha }}}^p$, in which a characteristic length appears to adimentionalize the curvature. By defining the free energy as such, the yield condition, which depends on the stress conjugate to the dislocation measure, depends on the curvature, therefore, the yield conditions become a partial differential equation in which there exists a dependency on the second gradient of the displacement field.

Regarding polycrystals, they adopted another measure to quantify the incompatibility of the plastic field:

$${\underset{\sim }{\mathsf{\eta }}}^p=\mathrm{curl}\left(\mathrm{curl}{\underset{\sim }{\mathsf{\epsilon }}}^p\right);$$

(39)

where $\underset{\sim }{\mathsf{\epsilon }}$ is the usual strain. Once again, the Clausius–Duhem inequality is explicitly given as summation of terms, and one of them depends on this incompatibility measure. By identifying the stress conjugate to this measure, and by giving the yield condition as function of this stress, it is possible to verify that the yield condition contains fourth order gradient of the displacement field.

Among the SGTs incorporating dislocation measures, the research line initiated by Wulfinghoff et al. also deserves to be mentioned. This viscoplastic theory was originally developed in [97], and then successively modified to incorporate the effect of grain boundary resistance on the plastic flow [98]. They adopted the combination of a micromorphic-like approach (see Section 3.9) with a penalty term in the free energy to incorporate the effect of the plastic strain and its gradient on the yield limit development.

3.8. Gradient of Micro-Structure Rotation—Cosserat Media

The Cosserat model derives from the work of the Cosserat brother at the beginning of the 20th century [99]. This continuum mechanics theory enhances the kinematics of the continuum with additional degrees of freedom, that are meant to represent the microstructural rotation. In the three-dimensional space, the Cosserat theory possesses three additional degrees of freedom: the three components of the pseudo-vector $\underline{\mathit{\omega }}$ that represents the local rotation of the microstructure. This vectorial field is defined, as the displacement field, over the entire domain under consideration. Another additional feature of the Cosserat kinematics is the definition of the gradient of the vector of microrotation, $\underline{\mathit{\omega }}\otimes \underline{\mathbf{\nabla }}$, which is sometimes referred to as wryness, and it incorporates both the curvature and the torsion of the microstructure. Based on the defined kinematic measures, the internal power reads:

$$p^{\left(i\right)}=\underset{\sim }{\mathsf{\sigma }}:\dot{\underset{\sim }{\mathbf{e}}}+\underset{\sim }{\mathsf{\xi }}:\dot{\underset{\sim }{\mathsf{\Gamma }}};$$

(40)

where $\underset{\sim }{\mathbf{e}}$ and $\underset{\sim }{\mathsf{\Gamma }}$ are the Cosserat strain and wryness, defined as:

$$\begin{array}{c}\underset{\sim }{\mathbf{e}}=\underline{\mathbf{u}}\otimes \underline{\mathbf{\nabla }}+\mathsf{\Omega };\end{array}$$

(41)

$$\begin{array}{c}\underset{\sim }{\mathsf{\Gamma }}=\underline{\mathit{\omega }}\otimes \underline{\mathbf{\nabla }};\end{array}$$

(42)

where $\underset{\sim }{\mathsf{\xi }}$ is the couple stress conjugate of the wryness, and the tensor $\underset{\sim }{\mathsf{\Omega }}$ is obtained from the micro-rotation vector $\underline{\mathit{\omega }}$ through the Levi-Civita permutation operator:

$$\underset{\sim }{\mathsf{\Omega }}=-\frac{1}{2}\underset{\simeq }{\mathsf{ϵ}}·\underline{\mathit{\omega }};$$

(43)

which returns a second order tensor. To be noted that if material and microstructure rotations are equivalent, the Cosserat strain coincides with the classical strain, the wryness becomes equal to the curvature in the Fleck1997 SGPT, and the Cosserat theory becomes a constrained Cosserat theory [100]. The Cosserat theory therefore employs the gradient of the rotation vector to regularize the boundary value problem. Assuming an elastic isotropic behavior, the constitutive model of the Cosserat media can be written as:

$$\begin{array}{c}\underset{\sim }{\mathsf{\sigma }}=\lambda \mathrm{tr}\left({\underset{\sim }{\mathbf{e}}}^e\right)\underset{\sim }{\mathbf{I}}+2\mu \mathrm{sym}\left({\underset{\sim }{\mathbf{e}}}^e\right)+2{\mu }_c\mathrm{skw}\left({\underset{\sim }{\mathbf{e}}}^e\right);\end{array}$$

(44)

$$\begin{array}{c}\underset{\sim }{\mathsf{\mu }}=\alpha \mathrm{tr}\left({\underset{\sim }{\mathsf{\Gamma }}}^e\right)\underset{\sim }{\mathbf{I}}+2\beta \mathrm{sym}\left({\underset{\sim }{\mathsf{\Gamma }}}^e\right)+2\gamma \mathrm{skw}\left({\underset{\sim }{\mathsf{\Gamma }}}^e\right);\end{array}$$

(45)

from which, two (more can be identified [64]) elastic characteristic lengths can be identified as [101]:

$$\begin{array}{c}{l}_{el}^b=\sqrt{\frac{\gamma }{\mu }};\end{array}$$

(46)

$$\begin{array}{c}{l}_{el}^t=\sqrt{\frac{\alpha +\beta +\gamma }{2\mu }};\end{array}$$

(47)

where the superscripts t and b indicate whether the characteristic length governs the torsional or bending response of the material respectively. The constrained Cosserat theory can also be obtained by enforcing the skew symmetric part of the Cosserat strain to be negligible through the employment of high enough values of the coefficient relating this term to the stress (${\mu }_c$) in Equation (44). The plastic framework of the Cosserat media can be taken from the one developed by de Borst [102], in which the J2 von Mises plasticity is extended to the Cosserat framework, and the equivalent stress measure reads:

$${\sigma }_{\mathrm{eq}}=\sqrt{\frac{3}{2}\left(\mathrm{a}{\underset{\sim }{\mathsf{\sigma }}}^{\prime }:{\underset{\sim }{\mathsf{\sigma }}}^{\prime }+\mathrm{a}{\underset{\sim }{\mathsf{\sigma }}}^{\prime }:{\underset{\sim }{\mathsf{\sigma }}}^{\prime \mathrm{T}}+\mathrm{b}\underset{\sim }{\mathsf{\mu }}:\underset{\sim }{\mathsf{\mu }}+\mathrm{b}\underset{\sim }{\mathsf{\mu }}:{\underset{\sim }{\mathsf{\mu }}}^{\mathrm{T}}\right)};$$

(48)

from which the plastic characteristic length of the Cosserat Media can be obtained as:

$$l_p=\sqrt{{\displaystyle \frac{\mathrm{a}}{\mathrm{b}}}};$$

(49)

3.9. Gradient of Micro-Structure Deformation—Micromorphic Media

The Cosserat medium allows the continuum to have an independent rotation from the material one. Similarly, the micro-strain framework allows the continuum to have an additional, independently defined strain (six additional degrees of freedom) that differs from the classical strain definition [78]. In reality, these two models are nothing but a constrained micromorphic model in which the micro-strain is constrained (Cosserat or Micropolar) or the micro-rotation is constrained (micro-strain). This theory, independently developed by Mindlin [52] and Eringen and Suhubi [103], adopts the gradient of the additional independent microdeformation field to describe the complete status of the domain, thus adding nine more independent variables to describe the media [104]. The internal power in the micromorphic media reads:

$$p^{\left(i\right)}=\underset{\sim }{\mathsf{\sigma }}:\left(\underline{\mathbf{u}}\otimes \underline{\mathbf{\nabla }}\right)+\underset{\sim }{\mathit{a}}:\underset{\sim }{\mathsf{\varphi }}+\underset{\simeq }{\mathbf{b}}⠇\left(\underset{\sim }{\mathsf{\varphi }}\otimes \underline{\mathbf{\nabla }}\right);$$

(50)

where $\underset{\sim }{\mathsf{\varphi }}$ is the second order micromorphic deformation tensor containing the additional nine degrees of freedom, $\underset{\sim }{\mathit{a}}$ and $\underset{\simeq }{\mathbf{b}}$ are power conjugates to the micromorphic deformation and to its gradient. In this framework, following Eringen [104], a single characteristic length arises naturally from the angular moment balance equation:

$$\underset{\sim }{\mathsf{\varphi }}-l^2{\mathbf{\nabla }}^2\underset{\sim }{\mathsf{\varphi }}=\underline{\mathbf{u}}\otimes \underline{\mathbf{\nabla }};$$

(51)

that links the displacements with the Micromorphic degrees of freedom [105].

The full micromorphic theory enhances the continuum with nine additional degrees of freedom, but it is still possible to take advantage of the versatility and the easiness of the numerical implementation of this theory by adopting a scalar modification of the micromorphic theory. Such solution was exploited by Anand et al., who used a single additional degree of freedom provided by the simplified micromorphic theory to model an equivalent measure of the cumulative equivalent plastic strain [106]. In this case, the micromorphic theory assumes a form similar to the previously discussed SGT, but it distinguishes itself from the others in the solid thermodynamic foundations on which it relies.

The micromorphic approach has also been used by Poh et al. as an implicit formulation that accommodates the gradient of the plastic strain [107]. The internal power that they defined has the following form:

$$p^{\left(i\right)}=\underset{\sim }{\mathsf{\sigma }}:\dot{\underset{\sim }{\mathsf{\epsilon }}}+\underset{\sim }{\mathbf{Q}}:{\dot{\underset{\sim }{\tilde{\mathsf{\epsilon }}}}}^p+\underset{\simeq }{\mathsf{\tau }}:\left({\dot{\underset{\sim }{\tilde{\mathsf{\epsilon }}}}}^p\otimes \underline{\mathbf{\nabla }}\right);$$

(52)

where $\underset{\sim }{\mathsf{\sigma }}$ and $\underset{\sim }{\mathsf{\epsilon }}$ are the usual Cauchy stress and strain, ${\underset{\sim }{\tilde{\mathsf{\epsilon }}}}^p$ is the microscopic plastic strain, $\underset{\sim }{\mathbf{Q}}$ and $\underset{\simeq }{\mathsf{\tau }}$ are the work conjugate to the microscopic plastic strain and its gradient. The microscopic plastic is a non-local measure of plastic deformation of the continuum and it is different form the the plastic strain. The governing equation of the microscopic plastic strain is obtained by developing the local equilibrium condition of the angular moment, that transforms into a governing equation of the higher degree of freedom, that is, the microscopic plastic strain, and it reads:

$${\underset{\sim }{\tilde{\mathsf{\epsilon }}}}^p={\mathsf{\epsilon }}^p+l_p^2{\underline{\mathbf{\nabla }}}^2{\underset{\sim }{\tilde{\mathsf{\epsilon }}}}^p;$$

(53)

where $l_p$ is a plastic characteristic length tuning the magnitude of the non-locality of the formulation. If $l_p=0$, the classical theory is retrieved. Equation (53) is a Helmholtz-type differential equation. The scalar version of this theory was previously developed by Peerlings [108], and it can be easily retrieved by adopting a scalar plastic microstrain instead of a tensorial one. In the scalar version, the definition of the free energy was given in a quadratic form as:

$$u^{\left(i\right)}={\displaystyle \frac{1}{2}}{\underset{\sim }{\mathsf{\epsilon }}}^e:\underset{\approx }{\mathbf{C}}:{\underset{\sim }{\mathsf{\epsilon }}}^e+{\displaystyle \frac{h}{n+1}}{p}^{n+1}+{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{3}{2}}\tilde{h}\right){(\tilde{p}-p)}^2+{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{3}{2}}\tilde{h}\right)l_p^2{\left(\mathbf{\nabla }\tilde{p}\right)}^2;$$

(54)

where $\underset{\approx }{\mathbf{C}}$ is the elastic fourth-order tensor, h and $\tilde{h}$ are hardening moduli, incorporated in the is free energy to model hardening, p and $\tilde{p}$ are equivalent plastic strain and equivalent plastic microstrain, n a material parameter.

4. Reported Applications of SGPTs in Manufacturing Processes Simulation

The implementation of a SGPT in a FEM solver is not an easy task, and the more complex the model, the more difficult is the implementation. Furthermore, if the SGPT has to be applied to machining simulation, this represents an additional element of difficulty due to the additional features that must be taken care of, such as friction, temperature and strain rate influence, large deformation framework and so on.

From the point of view of FE implementation, two groups of SGPT can be delineated: one that only requires to access and enhance the material behavior of the Finite Elements, and another one that requires a deeper level of accessibility. For the latter, a total new Element formulation must be provided. The key feature to differentiate between the two groups is the necessity of defining new degrees of freedom for the element. Whether it is demanded to add another degree of freedom or not, it depends on the nature of the continuum description that must be used: if the strain gradient needs to be evaluated through solution of PDEs, then the variable carrying the strain gradient measure requires to be treated as an additional degree of freedom, otherwise, the strain gradient will be treated as a solution-dependent variable. In the latter case, the strain gradient depends on the displacement field on a similar way the classical strain measure does, and the classical structure of the continuum description is unaltered.

Among the SGT previously presented, the CMSGT is the only one that does not require modifications of the element formulation. This traduces in a easier implementation of the theory, and this represent an advantage that should not be superficially neglected, especially if combined with the already mentioned difficulties characterizing machining simulations. Nonetheless, on a similar low-order theory developed by Acharya and Bassani [109] and Bassani [110], doubts have been raised on the physicality of the obtained response and on the well-posedness of the boundary value problem.

In the present section, several applications of SGPTs to machining operations will be presented. Although it is possible to find plenty of researches in literature focusing on the implementation and verification of SGPTs, here the objective is to gather the main contributions toward the application of SGPTs for manufacturing processes. Important advancements were done in the theoretical understanding of the phenomenon occurring during machining operations [5,6], but unfortunately not many contributions on the implementation of these SGPT for machining were found.

4.1. Scalar SGPT Applied to Flat Punch Molding

The first example of applications of SGPT in manufacturing operations is the one provided by Guha et al. [111]. They implemented the Fleck1997 SGPT as revisited by Fleck and Hutching in 2001 (ref. Section 3.4) in a viscoplastic isothermal large deformation framework to simulate flat punch molding. This manufacturing process is similar to the indentation test in case the size of the surface of the punch is comparable to the material grain size, thus achieving similar deformation fields to the ones developed during an indentation test. This test is of particular interest for manufacturing processes acting at macroscale ($\approx \mathsf{\mu }$m).

The main disadvantage of this SGPT is that the equivalent plastic strain is treated as an independent variable of the continuum, therefore it requires to be embedded as an additional degree of freedom in a Finite Element Procedure. Consequently, the independent field of the equivalent plastic strain must be solved through PDEs, and since the plastic deformation is not intrinsically defined over the domain (but rather it derives from the current deformation state of the media) the domain of definition of the equivalent plastic strain field requires to be updated at every iteration. The procedure of finding the boundaries of the definition domain of the equivalent plastic strain is facilitated if a Corotational approach or an Updated approach are used.

They performed simulations of a flat punch molding using different values of characteristic length. A strong dependency of the characteristic hardness on the characteristic length (or on the punch width) is predicted, as reported in Figure 9. As expected, this trend is very similar to the one characterizing the indentation tests [23,24,112]. Their research further emphasize the necessity of adopting a SGPT to simulate manufacturing operations acting on the micro scale.

4.2. SGPT Used to Model Rolling at Small Scale

The employment of the full gradient of plastic strain for manufacturing operations is reported in a recent contribution by Nielsen et al. [113]. They used the SGPT described by Gudmunson [81] first, and Fleck and Willis [114] later, to simulate rolling processes of metal sheets of up to 5 $\mathsf{\mu }$m of thickness in 2D. Simulations were performed under steady-state condition, integrations were done with forward Euler scheme, and the continuum was enhanced with nine additional degrees of freedom, being the components of the plastic strain tensor. Gaussian integration schemes were used for integration over the element. The contact between the plate and the roll was modeled through a sticking-sliding model. The material hardening was set to be dependent on the strain gradient.

Rolling simulations were carried out using different ratios of the characteristic length over the sheet thickness, and the results replicated the experimental values accurately. The simulated rolling force demonstrated sensitivity over the size-effect, similarly to what was experienced experimentally. The authors outlined that large size-effects were measured when the characteristic length over thickness ratio was 0.25.

4.3. MSGT Applied to Orthogonal Cutting

The first line of research here reported is the one originated from Liu et al. [8,115,116]. They implemented the Mix-Gao MSGT, thereby adopting a non-local Taylor-derived material hardening description as in Equation (34), where the equivalent measure of the strain gradient was evaluated as in Equation (33), and they used a non-local volume integral to evaluate the strain gradient at every integration point as proposed by Gao and Huang [117]:

$${\eta }_{ijk}={\int }_V\left[{\epsilon }_{ik}{\xi }_j+{\epsilon }_{jk}{\xi }_i-{\epsilon }_{ij}{\xi }_k-{\displaystyle \frac{1}{4}}\left({\delta }_{ik}{\xi }_j+{\delta }_{jk}{\xi }_i\right){\epsilon }_{pp}\right]dV{\left({\int }_V{\xi }_m{\xi }_{k}dV\right)}^{-1};$$

(55)

where ${\xi }_i$ is the i-th local coordinate over which integration is performed. In this framework, the classical continuum mechanics description stays the same, and the strain gradient effects only the plastic flow, acting as an additional source of hardening (this is after all the reasoning behind the MSTG). The flow stress they used is the one in Equation (34), and they used the Johnson–Cook evolution of the flow stress to replace f [118]. The characteristic length was identified with:

$$l=18{\alpha }^2{\left({\displaystyle \frac{\mu }{{\mathsf{\sigma }}_Y}}\right)}^{2}b;$$

(56)

They implemented this model in a Finite Element Procedure, using ABAQUS^® Implicit through a UMAT FORTRAN subroutine, and they used this model to study the effect of the radius of the cutting tool during orthogonal micro-cutting [8], and to verify the increments in material hardening and cutting energy by using the SGPT [116]. In Figure 10a a comparison is reported from the work of Liu et al. in terms of Specific Cutting energy while using different tool radii. In Figure 10b instead, the effect of using the SGPT is reported again in terms of Specific Cutting Energy. They used a sticking-sliding friction model as the one developed by Zorev [70]. The thermal fields during simulations have been solved through the classical heat equation. Heat was assumed to generate at the tool-workpiece contact location due to friction and at the plastic zones due to plastic work; the Taylor–Quinney constant was set to 0.9 (90% of the plastic work would transform into heat) [119]. An adaptive remeshing algorithm was used to obtain chip separation. The minimum mesh size they used measured 0.06 $\mathsf{\mu }$m in the cutting zone, with a chip thickness of 0.5 $\mathsf{\mu }$m.

The work of Liu et al. represents a big contribution toward the application of SGPTs for machining simulation, both in terms of progress toward the objective and in terms of proof of the fact that SGPT are necessary to achieve it. However, the disadvantages of using a MSGT are that the system of PDE describing the continuum behavior have not been rigorously derived and thermodynamical consistency of the model has not been proved, therefore, for future development, a more complete and rigorous SGPT would be ideally used, or a thermodynamically-consistent MSGT would have be provided.

Another contribution in simulating machining with SGPTs is the one provided by Asad, Mabrouki et al. [120,121,122]. They also used the MSGT as the one used by Liu et al., but they estimated the equivalent strain gradient measure based on geometrical consideration of the workpiece and the tool used for machining operation [121]. This has large implications on the limited applicability of this theory for a generic simulation. They also quantified the characteristic length as in Equation (56). Besides the evaluation of the strain gradient, they used the same flow stress description as the one used by Liu et al., but the chip formation was induced through a damage-enriched Johnson–Cook behavior law, where damage initiation and propagation were defined based on ductile Johnson–Cook fracture behavior [120]. They also used an ABAQUS^® Explicit VUMAT subroutine to implement the model. In Figure 11 the difference between using this variety of MSGT and the classical continuum mechanics description is reported for different cutting speeds.

4.4. Micromorphic Media Applied to Forming

The discussed Micromorphic theory is enhancing the medium by nine additional degrees of freedom, this being an equivalent deformation tensor of the microstructure defined at each point of the continuum. It has also been said that the Micropolar and the Microstrain theories are nothing but a constrained version of the Micromorphic theory. Due to the physics that Eringen described with the Micromorphic theory, this terminology associates to the introduction of additional degrees of freedom in the system that are connected with the microstructure. Through the years, many theories that enhance the continuum with additional degrees of freedom have been semantically linked with the Micromorphic theory, although the physics they described was not related to microstructure. However, in the literature, this became an established nomenclature for continuum mechanics theories enhancing the media with additional degrees of freedom.

The additional degree of freedom can be used for different purposes, it can be used for carrying an equivalent measure of the cumulative plastic strain [62,106,107,108] or for an equivalent damage variable [59,64]. In any case, the equation governing the angular momentum (that for classical continuum mechanics ensures the symmetry of the Cauchy stress tensor) transforms into a balance equation for the additional degrees of freedom and it is referred to as generalized balance equation. In case these additional variables are meant to replace already existing variables, the generalized balance equation transforms into a coupling equation between additional and already existing variables.

In 2017, Diamantopoulou et al. used the Micromorphic theory to enhance the continuum with one scalar additional degree of freedom (meant to represent the damage variable of the medium) and they simulated some metal forming processes with this theory [123]. The enhanced internal power of the continuum reads:

$$p^{\left(i\right)}=\underset{\sim }{\mathsf{\sigma }}:\left(\dot{\underset{\sim }{\mathbf{u}}}\otimes \underline{\mathbf{\nabla }}\right)+Y\dot{d}+\underline{\mathbf{Z}}·\underline{\mathbf{\nabla }}\dot{d};$$

(57)

where $\underset{\sim }{\mathsf{\sigma }}$ is the Cauchy stress tensor, d is the additional Micromorphic damage variable, Y and $\underline{\mathbf{Z}}$ are higher order stresses associated to the damage variable and its gradient respectively. The description of the theory continues with the considerations related to the thermodynamical admissibility of the theory: the definition of the Helmholtz free energy, the definition of the dissipation potential and the subsequent derivation of the material constitutive behavior and the elastoplastic flow rules from the previously defined potentials. In this medium description, the characteristic length appears in the generalized balance equation coupled with the gradient of the Micromorphic damage variable.

The theory was implemented in ABAQUS^® Explicit, through a VUEL subroutine. The introduction of an additional degree of freedom implies a modification of the system of PDEs to be solved, thereby compelling the modification of the Finite Element description rather than an higher-level modification of the material module.

The authors verified the developed formulation through tensile test and bending both performed using the material parameters of the DP1000 steel. The results of their analysis demonstrated an embedded mesh regularization property of the model, that was proven by varying the mesh size, whereas the global force-displacement graphs remained unchanged.

5. Conclusions

By distinguishing between higher order and higher grade formulations, the choice of implementing a continuum mechanics theory belonging to one group or the other depends on the specifications of the process that must be simulated. In general, higher grade theories are more complex to be implemented from a FEM point-of-view, whereas higher order theories are, in most cases, an extension of the classical FEM implementation to an higher number of degrees of freedom, thereby they result more straightforwardly implementable. In

Table 1. Summary of the presented SGPTs.

Name	# d.o.f.	# Charact. Lengths	Additional Feature	Higher Order/Higher Grade
Aifantis [36,79,80]	3	1/2	${\sigma }^Y\propto {\nabla }^2{\epsilon }_{eq}^p$	Higher Grade
Fleck1993 [85]	3	2	$p^{\left(i\right)}\propto \mathrm{skew}\left(\dot{\underline{\mathbf{u}}}\otimes \underline{\mathbf{\nabla }}\right)$	Higher Grade
Fleck1997 [20]	3	6	$p^{\left(i\right)}\propto \underline{\mathbf{u}}\otimes \underline{\mathbf{\nabla }}\otimes \underline{\mathbf{\nabla }}$	Higher Grade
Fleck2001 [42]	3	1	$p^{\left(i\right)}\propto {\mathsf{\epsilon }}^p\otimes \underline{\mathbf{\nabla }}$	Higher Grade
Gurtin [89]	3	2	$p^{\left(i\right)}\propto ({\underset{\sim }{\mathsf{\epsilon }}}^p\otimes \underline{\mathbf{\nabla }},\underset{\sim }{\mathsf{\theta }})$	Higher Grade
Gudmundson [81]	3	6	$p^{\left(i\right)}\propto {\mathsf{\epsilon }}^p\otimes \underline{\mathbf{\nabla }}$	Higher Grade
Gao [90,91]	3	1	${\sigma }^Y\propto \underline{\mathbf{u}}\otimes \underline{\mathbf{\nabla }}\otimes \underline{\mathbf{\nabla }}$	Higher Grade
Cosserat [99]	6	3	$p^{\left(i\right)}\propto (\underset{\sim }{\mathsf{\theta }},\mathsf{\Gamma })$	Higher Order
Micromorphic [103]	12	2	$p^{\left(i\right)}\propto (\underset{\sim }{\mathsf{\varphi }},\underset{\sim }{\mathsf{\varphi }}\otimes \underline{\nabla })$	Higher Order

the main characteristics of the presented SGPTs are reported.

For the simulation of machining operations that involve the development of high temperatures, regardless of the choice between higher order or higher degree, it is require the adoption of a theory incorporating a thermodynamical description that would ensures a proper thermodynamically-consistent description. The suggested approach would be to define the dissipation potential and the Helmholtz free energy from which the constitutive behavior of the medium would directly follow. If higher temperatures are foreseen, the dissipation potential assumes a fundamental role in determining the material response, thereby a proper calibration of the potential at operational conditions is required. If extremely high temperature develops, e.g., through the friction stir welding technique, melting temperature is likely to occur, and the solid description of the medium is not adequate anymore, since the continuum behaves as liquid rather than solid, therefore a mixed Lagrangian-Eulerian formulation is required to performed the simulation. However, the effect of the strain gradient during friction stir welding is not the main character of the whole process, therefore this manufacturing technique is left out of this manuscript.

One of the main drawback of the majority of the higher order theories is that, since the plastic strain is an additional variable of the system, the boundary conditions in terms of plastic strain must be applied at the moving elasto-plastic boundary, so another additional step is required to identify this boundary inside the original domain. Furthermore, if the simulated manufacturing technique involves the development of plastic strain at the domain boundaries, e.g., orthogonal cutting, then the choice related to the type of boundary conditions (in terms of plastic strain) that must be applied at this boundary is very important. If, however, the process that must be simulated does not involve the production of plastic strain at the boundary, then this question does not need to be posed, and the choice of using this type of continuum theory is not unfavorable. Two simulations of orthogonal cutting using SGPTs were reported, and both of them are using MSGT. In light of the previous discussion, the choice of the MSGT is a smart solution, since it has already been mentioned that this type of SGPT does not require the plastic strain field to be solved through PDEs, therefore it is not characterized by this type of problem.

To conclude, the choice of the SGPT to be adopted for simulations of manufacturing operations must be taken while looking at several characteristic of the observed phenomena occurring during the process:

• The forces/moments used to shear/shape the continuum localize in areas that are comparable to the grain size of the metal; in this case a further hardening due to the dislocations movement must be properly captured by adopting a SGPT;
• The velocity with which the continuum is deformed might induce material to behave viscously, and this must be covered by a proper design of the material behavior;
• The temperature developed during the process might be relatively high, and this imposes to considerate the thermal behavior (often adiabatic) that are taken into account by the definition of a thermodynamically-consistent description of the continuum;
• Implementation-wise if the plastic deformation is expected to occur at the boundaries, a higher degree theory is favorable.