The construction of regularization operators presented in this work is based on the introduction of strain or damage micromorphic degrees of freedom in addition to the displacement vector and of their gradients into the Helmholtz free energy function of the constitutive material model. The combination of a new balance equation for generalized stresses and of the micromorphic constitutive equations generates the regularization operator. Within the small strain framework, the choice of a quadratic potential w.r.t. the gradient term provides the widely used Helmholtz operator whose regularization properties are well known: smoothing of discontinuities at interfaces and boundary layers in hardening materials, and finite width localization bands in softening materials. The objective is to review and propose nonlinear extensions of micromorphic and strain/damage gradient models along two lines: the first one introducing nonlinear relations between generalized stresses and strains; the second one envisaging several classes of finite deformation model formulations. The generic approach is applicable to a large class of elastoviscoplastic and damage models including anisothermal and multiphysics coupling. Two standard procedures of extension of classical constitutive laws to large strains are combined with the micromorphic approach: additive split of some Lagrangian strain measure or choice of a local objective rotating frame. Three distinct operators are finally derived using the multiplicative decomposition of the deformation gradient. A new feature is that a free energy function depending solely on variables defined in the intermediate isoclinic configuration leads to the existence of additional kinematic hardening induced by the gradient of a scalar micromorphic variable.
Regularization differential or integral operators are widely used in mechanics in order to smooth discontinuities or restore the well-posedness of boundary-value problems. They are phenomenologically introduced in convection–diffusion problems of fluid mechanics and heat transfer  and in the damage of solids , with close links to filtering techniques in image analysis . Discontinuous variables like plastic strain in conventional plasticity can be smoothed at boundaries and interfaces to better reproduce physical deformation mechanisms, or in strain localization bands in the case of softening mechanical behaviour. Sharp interfaces are replaced by smooth interfaces in phase field models to simulate moving boundaries and thus avoid complex front-tracking methods [4–6]. The regularization is used in the two latter cases to obtain discretization-objective simulations results, i.e. fields that do not depend on the finite element, finite difference or Fourier grid size.
Various types of regularization methods are available relying on non-local integral operators , gradient formulations  or extra-degrees of freedom for smoothing strain or damage fields . The last two techniques very often involve so-called Helmholtz-type partial differential equations (PDEs) including a Laplace operator  or even bi-Helmholtz operators , see  for an anisotropic Helmholtz operator. Diffusion-like operators are used in phase field models in the form of Ginzburg–Landau or Allen–Cahn equations. Close relations exist between regularization operators used in continuum damage mechanics and in phase field theory, as recognized recently [12–15].
The regularization is often seen as a mathematical tool to be introduced into the physical model in an heuristic or ad hoc manner. It is then presented in the form of computational recipes to enhance existing algorithms. This is for instance the case of the regularization strategy proposed for softening plasticity and damage in [9,16–19] where the additional field equations and their coupling to the physical equations are postulated as PDEs independent of the specific form of the mechanical constitutive equations. By contrast, thermodynamically based formulations have been proposed where regularization differential operators are derived from new balance equations of generalized forces [20–23]. The form of the regularization operator then follows from the choice of constitutive equations linking generalized stresses and generalized strains and is not given a priori. In particular, the theories resorting to extra-degrees of freedom and of their gradients can be formulated within the micromorphic approach to gradient elastoviscoplasticity, as proposed in , the name micromorphic denoting additional strain-like degrees of freedom describing microstructure evolution according to Mindlin  and Eringen & Suhubi . The proposed procedure is especially useful to derive the form of the regularization operators under anisothermal or multiphysics conditions . The thermodynamic foundations of phase field models are well established  and a unified thermomechanical formulation is possible to reconcile the gradient, micromorphic and phase-field model classes [21,27–29].
The enhancement of the free energy density function by gradient terms of strain, damage or phase-field variables is usually limited to a quadratic contribution leading to a linear relationship between the generalized stresses and the gradient terms. Most applications are also limited to the small strain framework. They deal not only with problems of strain and damage localization in the mechanics of materials and structures, but also more recently with biomechanics and energy storage [30,31]. The linearized structure of these theories ensures in most cases good regularization properties of Helmholtz or diffusion operators. The formulation and performance of such operators in nonlinear cases like non-quadratic dependence or large deformations remain largely open.
First, extensions of regularization approaches to finite deformations dealt with the simulation of strain localization in gradient materials [32,33] or non-local media  where Lagrangian and Eulerian formulations were presented and compared. The approach based on extra-degrees of freedom was first extended using an Eulerian formulation, a logarithmic elastic strain and an expression of the Helmholtz equation involving the Laplace operator w.r.t. to the current configuration . Generalized Helmholtz equations governing micromorphic-like variables were proposed within the Lagrangian framework in [36,37] with quadratic potentials w.r.t. the Lagrangian gradient of the extra-degrees of freedom. A constitutive formulation using local rotating frames, i.e. hypoelastic laws, was developed in [17,38]. The authors in  made used of a multiplicative decomposition of the deformation gradient.
Guidelines for the design of regularization operators at large deformations can be deduced from a hierarchy of micromorphic and strain gradient constitutive models as proposed in [40–45]. In particular, decompositions of higher order deformation measures into elastic and plastic parts were proposed. These constitutive settings are more general than the ones targeted in this work. They aim at the formulation of physically based higher order constitutive equations and not only at regularization purposes. The aim of the present class of models is different, it concentrates on regularization operators and therefore focuses on extensions of standard elastoviscoplasticity models involving strain-like or hardening/softening variable-like additional degrees of freedom. The enhanced constitutive equations should be kept as simple as possible. In particular, dissipative contributions of higher order variables are not considered as a first step as such simplified models already provide powerful regularization properties. More specific dissipative higher order constitutive equations can be found in [40,41,46] and in the references quoted therein.
Motivations for nonlinear potentials with respect to gradient of the extra-degrees of freedom stem from recent strain gradient crystal plasticity models making use of rank one, power law or logarithmic functions of the gradient terms for better description of dislocation behaviour [47–50]. Non-quadratic potentials are seldom in the phase field community. Examples can be found for the modelling of grain boundary migration and recrystallization [51,52]. Such nonlinear relations between generalized stresses and gradient terms lead to fundamentally new kinds of differential operators whose properties are essentially unknown.
The objective of this article is to review and propose nonlinear extensions of previous micromorphic and strain/damage gradient models along two distinct lines: the first one introducing nonlinear material relations between higher order stresses and strains; the second one envisaging different classes of finite deformation formulations of the models. Following both directions, nonlinear regularization operators, with up to now mostly unknown mathematical properties, will be exhibited. We insist here on the generic character of the proposed approach, applicable in a systematic way to a large class of elastoviscoplastic and damage models allowing for anisothermal and multiphysics coupling.
The original micromorphic approach used to generate Helmholtz-like regularization operators is first recalled in §2 in the case of scalar and tensor microstrain degrees of freedom. Nonlinear extensions of the relation between generalized stresses and strains are proposed in §3, still within the small strain hypothesis. Three methods of extension to finite deformations are presented then: a fully Lagrangian method in §§4 and 5, the use of local objective rotating frames in §6 and finally in combination with the multiplicative decomposition of the deformation gradient (§7) which represents the best-suited method for anisotropic materials. In the final conclusions, the new contributions of this article are pointed out, together with the limitations and need for future developments of the approach.
An intrinsic notation is used throughout where zeroth, first, second, third and fourth order tensors are denoted by and , respectively. The simple, double and triple contractions are written ⋅, : and ⋮, respectively. The proposed notation can be translated into a direct index one with respect to an orthonormal Cartesian basis . For example, where repeated indices are summed. Tensor product is denoted by ⊗ and the nabla operator is ∇. For example, the component ijk of is aij,k, the comma denoting partial differentiation with respect to the ith coordinate. In particular, Δ=∇⋅∇ is the Laplace operator.
2. The original approach based on quadratic potentials
The micromorphic approach delivering linear Helmholtz-type regularization operators is illustrated for three kinds of micromorphic degrees of freedom, namely a microstrain tensor in the spirit of Eringen's theory, and scalar variables associated with equivalent total or plastic strain measures, respectively.
(a) Microstrain tensor model
According to Eringen's and Mindlin's original approach [24,25], the material point is endowed with the usual translational degrees of freedom (d.f.), the displacement vector, , and an additional microdeformation tensor describing the rotation and distortion of a triad of directors attached to the microstructure: . The theory is presented here in the case of a symmetric microstrain tensor following [20,41,53] which is less general than Eringen's original non-symmetric microdeformation. The microstrain components χij are generally distinct from those of the classical small strain tensor εij=(ui,j+uj,i)/2. The method of virtual power is used to derive the balance equations of the theory, following . The power density of internal forces is assumed to be the following linear form: 2.1 in the small strain context. The Cauchy stress tensor is found to be divergence free in the absence of body forces. The generalized stress tensors and fulfil the following static balance equation: 2.2 with the associated Neumann boundary conditions , where the double force distribution is prescribed at the boundary of normal . The dual Dirichlet conditions amount to impose the microdeformation at the boundary. The boundary conditions provided in this work are not heuristic, they are derived from the method of virtual power following .
The Helmholtz free energy volume density is a function of the strain tensor, possible internal variables, α, the microdeformation tensor and its gradient.
The entropy principle of thermodynamics is adopted in its local form: .
Substituting the dependency of the free energy function on the chosen state variables leads to the Clausius–Duhem inequality: 2.3 For regularization purposes, it is enough to consider dissipation-free contributions of the higher order stress, thus adopting the following state laws: 2.4 the first relation being the classical hyperelastic law. The free energy density is split into two main contributions: 2.5 where ψref refers to any classical reference mechanical model and ψχ to the generalized contributions including the microdeformation gradient and a constitutive variable , function of strain, microstrain and/or internal variables, to be specified. An efficient model is obtained by adopting a quadratic potential ψχ: 2.6 where fourth and sixth rank tensors of higher order moduli have been introduced. The constitutive choice, , leads to a coupling of strain and microstrain tensors in the constitutive model. The state laws (2.4) are then linear stress–strain relations: 2.7 The classical Hooke tensor has been introduced highlighting the coupling between strain and microstrain, being the plastic strain tensor in the decomposition of the total strain into elastic and plastic parts. The regularization operator is obtained after substitution of the state laws into the extra-balance equation (2.2): 2.8 It can be made more apparent when the simplifications are adopted, thus reducing the number of higher order moduli to two, Hχ and A. Then, the PDE (2.8) reduces to 2.9 where Δ is the Laplace operator. The characteristic length, , arises in the theory. As a result, the microstrain is solution of the PDE (2.9) of Helmholtz type where the classical strain tensor acts as a source term. Convexity of the free energy density function is assumed, thus ensuring the stability of the model. This implies the positivity of the higher order moduli, Hχ>0, A>0. This is an essential property of the model ensuring regularization properties.
In the simplified version, the coupling between strain and microstrain is apparent in the generalized Hooke law: 2.10 Very high values of the coupling modulus Hχ enforce the internal constraint , so that 2.11 The micromorphic model then reduces to a gradient elasticity model closely related to Aifantis so-called gradelas model [55,56].
(b) Scalar microstrain
A variant of the previous model consists in reducing the number of extra-degrees of freedom from 6 to 1: 2.12 where χ is a scalar microstrain variable. In the microdilatation model described in [41,57,58], this degree of freedom is related to the trace of Eringen's microdeformation tensor accounting for instance for microvoid volume changes. The generalized stresses of the theory then reduce to the scalar a and the vector related by the following extra-balance equation, superseding (2.2): 2.13 with the associated scalar Neumann boundary condition and the dual Dirichlet scalar condition prescribing χ at the boundary.
The two last state laws in (2.4) are replaced by the lower order ones: 2.14 The proposed more specific form of the free energy potential is 2.15 where the relative strain measure e:=εeq−χ involves the equivalent strain measure, εeq, function of the three invariants of the strain tensor . A quadratic form is 2.16 where ψα related to hardening is left unspecified and can be non-quadratic. Accordingly, 2.17 so that the scalar microstrain is solution of the following PDE that results from (2.13) and (2.17): 2.18 In the isotropic case , the regularization operator then takes the form 2.19 This operator was introduced as early as  using such a scalar microstrain variable related to an equivalent strain measure. The coupling between macro- and microstrain leads to an additional contribution to the stress tensor in the form 2.20 If, for instance, the Euclidean norm of the stress tensor is introduced as an equivalent strain measure, 2.21 the coupling term in (2.20) is nonlinear.
(c) Scalar plastic microstrain
As an alternative, the following relative strain measure is adopted in equation (2.15): e:=p−χ, where p is the cumulative plastic strain arising in the plasticity theory in such a way that the plastic power reads: . The equivalent stress σeq is the stress measure involved in the yield function: 2.22 The cumulative plastic strain p is used here as an internal variable, together with , to model isotropic hardening in plasticity. The initial yield stress in tension is R0, from (2.22).
The choice of the following partly quadratic potential 2.23 leads to the following state laws: 2.24 where Rref(p)=ρ∂ψα/∂p is the standard hardening law. It is apparent that the classical Hooke law is unaffected, in contrast with (2.20), whereas the hardening law is modified by the coupling between plastic macro- and microstrain.
It was noted in  that the regularized hardening law (2.24) differs from the function R(χ) which is obtained by substituting the variable p by the micromorphic/regularized variable χ in the hardening function, as initially proposed in . This inconsistency leads to difficulties in combining hardening and softening material behaviour, as recognized in [20,60,61], difficulties that are settled by adopting the previous thermodynamic framework.
The hardening rule can be related to higher order space derivatives of the plastic microstrain by combining the state laws (2.24) with the balance equation (2.13): 2.25 the latter expression being written for the isotropic case. This PDE is associated with Dirichlet boundary conditions prescribing the microstrain χ or Neumann conditions giving the flux at the boundary. Note that the hardening law is valid for hardening or softening materials depending on the sign of the plastic tangent moduli ∂Rref/∂p. When the internal constraint e≡0⇔p=χ is enforced, or, equivalently when the penalty modulus Hχ is high enough, the Laplacian of the cumulative plastic strain itself appears in the hardening rule, which corresponds to Aifantis celebrated strain gradient plasticity model [62,56]. Strain gradient plasticity there arises as a limit case of the micromorphic model. Note that the sign of the material parameter A in (2.25) was heavily debated in heuristically introduced higher order terms in Aifantis-like models. The present thermodynamic approach requires A>0 to ensure the convexity of the gradient term in the free energy, and therefore, the model regularizing power.
When the material point is under plastic loading conditions, the combination of the yield function (2.22) and of the hardening law (2.25) provides the current value of the equivalent stress measure 2.26 Let us give a simple one-dimensional example where σeq is homogeneous. Considering the limit case e≡0 and a simple linear hardening rule, Rref(p)=Hp, H being the plastic modulus, the cumulative plastic strain is solution of the differential equation 2.27 which delivers two types of solutions depending on the sign of the plastic modulus.
— Hardening materials, H≥0: the plastic field is of hyperbolic/exponential type plus a possible parabolic contribution coming from the constant term. This model describes size dependent boundary layer effects related to size effects in the behaviour of metals for instance. The found size effects are discussed in .
— Softening materials, H<0: the plastic field is of harmonic type in addition to the possible parabolic contribution. This corresponds to the localization of plastic strain into a band of finite width in the form of a sinus arch. This regularization property has been widely used, for instance, in the simulation of strain localization phenomena like shear banding [64,39].
Both types of solutions arise in the description of propagating localization bands as encountered in the Lüders phenomenon in steels [65,66]. They illustrate the regularizing character of the original model: smoothing of discontinuities at interfaces and boundary layers in hardening materials, and description of finite width localization bands in softening materials. The proposed micromorphic approach in the present isotropic version introduces only two additional parameters. Parameter A can be identified from the scaling laws for size effects in hardening plasticity or from the characteristic width of strain localization bands for softening behaviour. The parameter Hχ can be seen as a penalty term high values of which lead to the underlying gradient model, or as a true micromorphic constitutive parameter that can be identified again to better describe the size-dependent structural responses, as done in  for micromorphic crystal plasticity. The anisotropic case involving a large number of additional parameters is much more challenging from the identification perspective. Extended homogenization procedures can be used to identify the whole set of parameters from the consideration of an underlying periodic microstructure, as reviewed in [68–70].
3. Nonlinear strain gradient potentials
Recent results in the plasticity of metals have revealed the limitations of a quadratic potential ψ∇ in equation (2.15), especially regarding the scaling of size effects in dislocation plasticity [47–50,63,71,72]. Two classes of nonlinear potentials are explored below motivated by the latter references. The motivations come from crystal plasticity theory mainly but the models are developed here in the context of phenomenological laws for polycrystalline materials. In this section, the specific form of the studied free energy potential is 3.1 where χ is a plastic microstrain degree of freedom. A normalized anisotropic norm of its gradient is introduced in Appendix A: 3.2 involving the characteristic length ℓc and a definite positive symmetric second-order tensor (dimensionless).
(a) Power law potential
The following power law gradient potential is investigated: 3.3 where the power m≥1 for reasons of convexity. The elastic shear modulus μ sets the physical dimension of this energy contribution and ℓc is the actual constitutive length of the model, although other lengths could be defined from the ratio between and any combination of elastic or plastic moduli. The higher order stress is computed from the hyperelastic law (2.14) 3.4 where the direction is given in Appendix A by equation (A.2). Its divergence enters the balance equation (2.13) and is computed as 3.5 The plastic microstrain is then governed by the following PDE defining the differential operator Op: 3.6 and 3.7 The operator is anisotropic due to the tensor and nonlinear because of the presence of . The nonlinearity disappears in the case of a quadratic potential ψ∇, m=2, as in §2c. The original regularization operator (2.19) is retrieved in the isotropic case for the power m=2: 3.8 The case m=1 is of particular importance as it has been considered at various places in the literature [47,48,50]: 3.9 The operator is singular at g=0 due to the non-differentiability of the potential ψ∇, and additional regularization techniques must be used in order to perform finite-element simulations in this context [50,73]. The coupling between plasticity and the micromorphic variables takes place at the level of the hardening law, with the hardening variable α:=p: 3.10 Taking equation (3.5) into account provides the enhanced hardening law: 3.11 In the isotropic case, the previous expression specializes to 3.12 where the direction is the direction of the gradient of the scalar microstrain, according to equation (A.1). In the previous hardening law (3.12), the usual Laplace term is complemented by a non-linear contribution that vanishes for m=2. The special case m=1 for a rank one potential ψ∇ leads to the following hardening law deduced from (3.12): 3.13 The norm of the gradient of plastic cumulative plastic strain g can be regarded as the phenomenological counterpart of the norm of the dislocation density tensor in the physical crystal plasticity theory. The quadratic case m=2 was used in strain gradient plasticity as a first constitutive proposal [74,75]. It turns out that the size effects predicted based on a quadratic potential are not consistent with results from physical metallurgy [67,63]. The singular case m=1 provides consistent scaling laws for the yield stress as a function of channel width in laminate microstructures [47,50]. In particular, the singular character of the model results in a size-dependent abrupt increase of the apparent yield stress. By contrast, regular potentials lead to a size-dependent apparent hardening modulus. The reader is referred to  for applications of the power law model to size effects in hardening crystal plasticity.
(b) Logarithmic potential
Motivated by energy considerations in dislocation theory, several authors have considered logarithmic functions of scalar dislocation densities or of the norm of the dislocation density tensor [49,50,71,72]. In the present context of phenomenological metal plasticity, a logarithmic function of the norm of the gradient of the scalar plastic microstrain is proposed: 3.14 where , g still given by equation (3.2) and g0 is a constant. This function is convex with respect to g for g≥0, with ψ∇(0)=0. It is not differentiable at g=0. A possible regularization is to consider that the initial value of g is non-zero, for instance equal to g0. Other regularizing choices are possible like the use of a quadratic potential in the interval 0≤g≤g0 . The generalized stress vector 3.15 is directed along vector , its norm grows as a logarithmic function of g and vanishes at g=g0. The divergence of the generalized stress vector follows: 3.16 which leads to the following nonlinear operator: 3.17 In the isotropic case, it becomes 3.18 The enhanced hardening law takes the general form: 3.19 These nonlinear PDEs are associated with the same Neumann or Dirichlet boundary conditions as presented in §2c. The physical motivations for the prescription of such boundary conditions in the context of dislocation plasticity are related to the role of interfaces and passivated free surfaces as obstacles to dislocation motion (so-called micro-hard Dirichlet conditions). They have been discussed for instance in [76,77].
(c) One-dimensional example
In the one-dimensional case with a single non-vanishing stress component σ(x) (e.g. simple tension or shear) and with all variables depending on x only, the plastic microstrain variable is solution of the following differential equation derived from equation (3.7) for a power-law potential: 3.20 where the prime indicates derivation with respect to the single spatial coordinate. The enhanced hardening rule becomes 3.21 In the quadratic case, m=2, the linear regularization Laplace operator (2.19) and Aifantis-like model (2.25) are retrieved.
The case m=1 leads to the condition p=χ and no extra-hardening. This leaves the possibility of localization of plastic strain and plastic strain gradient in the form of interface dislocations as discussed in [49,50]. The corresponding singular distribution of plastic microstrain in conjunction with the relation (3.15) was shown in the latter reference to lead to a size-dependent overall increase of the apparent yield stress and to no extra-hardening.
The logarithmic potential (3.14), the associated regularization operator (3.18) and the enhanced hardening law (3.19) are now specialized to the one-dimensional case: 3.22 In the absence of classical hardening, i.e. when Rref(p)=R0, a constant initial threshold, the differential equation, χ′′=μℓc|χ′|, admits exponential solutions with cusps as illustrated in the example given in gradient plasticity in reference . The logarithmic potential is inspired from strain gradient crystal plasticity models emerging from the statistical theory of dislocations . It has the advantage with respect to the m=1 case that it can account for both enhanced strength and hardening at small sizes, as demonstrated in . It is expected that similar behaviour can be obtained from the present isotropic polycrystal model, as suggested by Ohno .
4. Micromorphic and gradient hyperelasticity
Nonlinearity arises not only from nonlinear material response but also from the consideration of finite deformations. The impact of finite strains on regularization operators is largely unexplored. It is first illustrated in the pure hyperelastic case, leaving aside for a moment the inclusion of plastic effects. The Lagrangian coordinates of the material points are denoted by on the reference configuration Ω0, whereas their positions in the current configuration Ω are called . The gradient operators with respect to Lagrangian (reference) and Eulerian (current) coordinates are denoted by ∇0 and ∇, respectively. The deformation gradient is , where the displacement vector is the function . The initial and current mass density functions are and , respectively.
(a) Finite microstrain tensor model
The linear microstrain model discussed in §2a is now extended to the finite deformation case by applying the general approach presented in . The additional degrees of freedom are the six components of a microstrain tensor , a second-order symmetric tensor associated with the right Cauchy–Green strain measure of the micromorphic deformation. Accordingly, the microstrain is regarded here as a Lagrangian variable. The Lagrangian power density of internal forces takes the form 4.1 p(i) being the Eulerian internal power density and the Jacobian . The right Cauchy–Green tensor is and is the Piola stress tensor. The method of virtual power can be used, see , to show that the generalized Lagrangian stress tensors, a second and a third rank tensor, fulfil the following balance equation: 4.2 in addition to the balance of momentum equation 4.3 in the absence of body and inertial forces, the divergence operator Div being computed with respect to Lagrangian coordinates. The corresponding Eulerian forms of the balance equations are 4.4 with the Neumann boundary conditions for tractions and double tractions: , involving the normal vector of the current surface element. The relations between the Lagrangian and Eulerian generalized stress tensors are 4.5 The hyperelastic free energy density function is and the stress–strain relations read 4.6 The relative strain tensor was defined in  and measures the difference between macro- and microstrain.
As an example and straightforward generalization of equation (2.6), the following potential is proposed: 4.7 where a single penalty modulus Hχ was introduced and where ψref is a standard hyperelastic strain energy potential (neo-Hookean, etc.). The higher order term involves a sixth-rank tensor of elasticity moduli which is symmetric and assumed definite positive . The stress–strain relations (4.6) become 4.8 Note that the classical hyperelastic relation is complemented by a coupling term involving the microstrain tensor. Taking the balance equation (4.3) into account, the set of PDEs for the microstrain components is found to be 4.9 The associated regularization operator is now given in the simplified case where the sixth rank tensor of higher order moduli is assumed to be the identity multiplied by the single modulus A: 4.10 It involves the Lagrangian Laplace operator Δ0(•)=(•),KK in a Cartesian frame where capital indices refer to Lagrange coordinates and the comma to partial derivation. It is linear w.r.t. Lagrangian coordinates but the full problem is of course highly nonlinear. The associated Eulerian partial differential operator is nonlinear in the form: χIJ,KK=χIJ,klFkKFlK, where small index letters refer to the current Cartesian coordinate.
An Eulerian formulation of the proposed constitutive equations is possible. It will be illustrated in the next section in the case of a scalar microstrain variable for the sake of brevity.
(b) Equivalent microstrain model at finite deformation
The set of degrees of freedom of the proposed model is given by equation (2.12) and contains a scalar microstrain variable . The latter variable is assumed to be a Lagrangian quantity, invariant w.r.t. change of observer. The free energy density is a function of the following argument . The corresponding state laws fulfilling the Clausius–Duhem inequality take the form 4.11 As an example and straightforward generalization of equation (2.6), the following Lagrangian potential is proposed: 4.12 The microstrain is compared with some equivalent strain measure, Ceq, function of the invariants of . The stress–strain relations (4.11) become 4.13 Taking the balance equation into account, the PDE governing χ is 4.14 Let us mention the corresponding isotropic form, when , 4.15 where Δ0 is the Laplace operator with respect to Lagrangian coordinates.
An example of equivalent strain measure which the microstrain is compared with is 4.16 thus generalizing the model (2.21) to finite strains.
The formulation of a constitutive law based on Eulerian strain measures, and ∇χ, is now envisaged. It relies on the choice of a free energy potential . Galilean invariance of the constitutive law requires that this function fulfils the following conditions: 4.17 for all constant orthogonal transformations . This amounts to stating isotropy of the function ψ. Representation theorems are available for such functions, , where the Bi are the eigenvalues of . The Cauchy stress tensor, , is known then to commute with such that 4.18 where the strain rate tensor, , is the symmetric part of the velocity gradient, .
The hyperelastic state laws then take the form 4.19 As an example and straightforward generalization of equation (2.6), the following Eulerian potential is proposed: where ψref refers to a standard isotropic elasticity potential in classical mechanics. Note that Beq=Ceq as and share the same eigenvalues. The state laws (4.19) become 4.20 The Eulerian regularization operator follows from (4.4): 4.21 where Δ is the Laplace operator with respect to the Eulerian coordinates.
It is essential to note that the isotropic regularization operators (4.15) and (4.21) are distinct. For, if the Lagrangian higher order elastic law is linear with respect to the constitutive quantities involved, the deduced Eulerian law is NOT linear: 4.22 so that the Eulerian regularization operator will not be linear, i.e. different from (4.21).
5. Finite deformation micromorphic elastoviscoplasticity using an additive decomposition of a Lagrangian strain
The most straightforward extension of the previous framework to viscoplasticity is to introduce a finite plastic strain measure in the decomposition of a Lagrangian total strain tensor. Such Lagrangian formulations of elastoviscoplasticity involve the additive decomposition of some Lagrangian total strain measure into elastic and viscoplastic parts: 5.1 Many choices are possible for the invertible function h with restrictions ensuring that is a strain measure (symmetric, vanishing for rigid body motions, differentiable at 0 so that the tangent is the usual small strain tensor , used before, see ). Hill's strain measures are obtained for power law functions such that , the latter being the Lagrangian logarithmic strain. The case m=2 corresponds to the Green–Lagrange strain measure for which this finite deformation theory was first formulated by Green and Naghdi, see [79,80] for the pros and the cons of various such formulations. This Lagrangian formulation is preferable to Eulerian ones based on corresponding Eulerian strain measures in order not to limit the approach to isotropic material behaviour . The additive decomposition of the Lagrangian logarithmic strain is put forward in the computational plasticity strategies developed in [82,83]. However, there is generally no physical motivation for the selection of one or another Lagrangian strain measure within this framework. This approach favours one particular reference configuration for which the corresponding strain is decomposed into elastic and plastic parts, again without clear physical argument. Changes of reference configuration lead to complex hardly interpretable transformation rules for the plastic strain variables, see  for a comparison of finite deformation constitutive laws with respect to this issue. Note also that limitations arise from using a symmetric plastic strain variable , especially when plastic spin relations are needed for anisotropic materials . This framework was applied to micromorphic and gradient plasticity and damage theories in [35,85,86].
A Lagrangian conjugate stress tensor is defined for each strain measure such that 5.2 The power density of internal forces is: , and the free energy density function has the following arguments: . The dissipation inequality then reads 5.3 from which the following state laws are selected 5.4 The flow and hardening rules can be determined from the suitable choice of dissipation potential : 5.5 The existence of such a dissipation potential is not necessary but assumed here for convenience. Alternative methods of exploitation of the dissipation principle exist for micromorphic continua, for instance, based on the extended Liu procedure [87,88].
Two straightforward extensions of the micromorphic approach to finite strain viscoplasticity based on an additive decomposition of a Lagrangian strain measure are proposed 5.6 where Eheq is an equivalent total strain measure, or, alternatively, 5.7 where is the cumulative plastic strain in the present context. These choices, respectively, provide the following regularization PDE: 5.8 If , then the regularization operator involves the Lagrangian Laplace operator Δ0 in the same way as in equation (4.10).
6. Finite deformation micromorphic viscoplasticity using local objective frames
An alternative and frequently used method to formulate anisotropic elastoviscoplastic constitutive equations at finite deformations that identically fulfil the condition of Euclidean form invariance (also called material frame indifference ), is to resort to local objective rotating frames, as initially proposed by Dogui & Sidoroff [89,90]. A local objective rotating frame is defined by the rotation field , objective w.r.t. to further change of observer, and taking different values at different material points and different times: 6.1 It is based on the idea that there exists for each material point a privileged observer w.r.t. which the constitutive law takes a simple form. The method is described in details in  and is used in many commercial finite-element codes with the standard choices: co-rotational frame, such that , being the skew-symmetric part of the velocity gradient, and polar frame, such that , being the rotation part in the polar decomposition of the deformation gradient . The main drawback of this method is the absence of thermodynamic background since, depending on specific constitutive choices within this framework, a free energy function of the strain may not exist.
This method is now applied to a micromorphic model including a scalar additional d.f. χ as in equation (2.12). Extension to higher order tensor-valued additional degrees of freedom is straightforward. The field equations are still given by (7.6). The stresses w.r.t. the local objective frames are 6.2 Time-derivation of these relations shows that the rotated stress derivatives are given by 6.3 i.e. objective derivatives of the corresponding Eulerian stress tensors. If the co-rotational frame is used, the corresponding time derivative is the Jaumann rate, whereas it is the Green–Naghdi stress rate when the polar rotation is used. The same procedure is applied to the strain rates 6.4 The time integration of the second equation in the rotating frame provides the variable . It must be underlined that is NOT equal to . It is NOT the exact material time derivative of a constitutive variable, in general. The standard procedure then consists in postulating an additive decomposition of the rotated strain rate into elastic and plastic parts as 6.5 where the elastic and plastic strain and , respectively, are solely defined in the rotated frame. Anisotropic elastic laws are assumed to take the form 6.6 Time-derivation of these equations and consideration of equation (6.2) show that the elasticity laws are in fact hypoelastic and that, generally, there does not exist a free energy density function from which they can be derived .
The yield function and the flow rule are formulated within the rotated frame 6.7 where normality is assumed for convenience and the viscoplastic multiplier is given by some viscoplastic law. The evolution of internal variables is of the form for suitable functions H. The yield radius is chosen as the following expression inspired by the previous thermodynamically based models: .
This extension of the micromorphic approach to finite deformations using rotating frames has been proposed first by  and used by these authors for metal forming simulations involving regularized damage laws. As a result, the regularization operator can be written as 6.8 In the isotropic case, , the regularization operator reduces to 6.9 It is worth insisting on the fact that, in general, . Accordingly, the previous equation does not involve the Laplace operator and the regularization is therefore nonlinear even with respect to rotated quantities.
Among all choices of rotating frames, the one associated with the logarithmic spin rate tensor  was claimed to be the only one such that, when , the isotropic hypoelastic strain–stress relation turns out to be hyperelastic. However, this property does not pertain to the general case , so that this specific choice does not in general provide any explicit form of the regularization operator.
Alternative constitutive choices are possible for the higher order stresses if Laplacian operators are preferred. They amount to restricting the use of the rotating frame only for the classical elastoviscoplasticity equations and to writing independently, , so that the regularization operator is expressed in terms of the Eulerian Laplace operator Δ, see (7.22) in the next section, or which leads to the Lagrangian Laplace operator Δ0, see (7.17) in the next section.
Note that limitations in the formulation of anisotropic plasticity arise from using symmetric plastic strain variable and that generalizations are needed in order to introduce necessary plastic spins for materials with microtructures, which is possible within the rotating frame approach .
7. Finite deformation micromorphic elastoviscoplasticity based on the multiplicative decomposition
The most appropriate thermodynamically based framework for the formulation of finite deformation anisotropic elastoviscoplasticity relies on the multiplicative decomposition of the deformation gradient, as settled by Mandel . This method is applied here to a generalized continuum model again limited to one scalar degree of freedom, χ, in addition to the displacement vector, , see (2.12). The gradients of the degrees of freedom can be computed with respect to the reference or current coordinates: 7.1 and 7.2 The consideration of microdeformation degrees of freedom of higher order is possible without fundamental modification of the approach below, see .
A multiplicative decomposition is envisaged in this section in the form 7.3 which assumes the existence of a triad of directors attached to the material point in order to unambiguously define the isoclinic intermediate local configuration, labelled (♯) in the sequel, see [94,95]. The directors are usually related to non-material microstructure directions like lattice directions in single crystals or fibre directions in composites. The existence of such directors is required for the formulation of objective anisotropic constitutive equations . The Jacobians of all contributions in equation (7.3) are denoted by 7.4 They are used to relate the mass densities with respect to the three local configurations.
In the present work, the microdeformation gradient is not split into elastic and plastic contributions, although it is possible as done in , at the expense of additional evolution laws to be determined and to drastically different regularization operators.
The power density of internal and contact forces are 7.5 The invariance of p(i) with respect to any change of observer requires the Cauchy stress to be symmetric. The scalar microstrain is assumed to be invariant. The application of the method of virtual power based on the previous densities leads to the following static field and boundary equations, in the absence of body forces: 7.6
(a) Lagrangian formulation
The Lagrangian free energy density is a function , where is the elastic strain and α a set of internal variables accounting for material hardening. Note that the usual elastic strain tensor is defined with respect to the intermediate configuration to comply with standard anisotropic plasticity, whereas is Lagrangian. The presented formulation is therefore not purely Lagrangian but rather mixed. The local Lagrangian form of the entropy inequality is . Accounting for the multiplicative decomposition (7.3), the power of internal forces is expanded as 7.7 where the Piola stress tensor w.r.t. the intermediate configuration and the Mandel stress tensor  are, respectively, defined as 7.8 The Lagrangian generalized stresses in (7.7) are a0=Ja and .
As a result the dissipation rate becomes 7.9 The following state and evolution laws ensure the positivity of : 7.10 and 7.11 Following Mandel , a dissipation potential , function of the driving forces for plasticity, is introduced to formulate the flow and hardening variable evolution rules. If the dissipation function is convex w.r.t. and concave w.r.t. R, the positivity of dissipation is ensured. Specific expressions for Ω within the context of viscoplasticity can be found in .
As an example, the following free energy potential is proposed: 7.12 where ψα is the appropriate free energy contribution associated with usual work-hardening (not specified here) and is the Green–Lagrange elastic strain measure. The microstrain variable is compared to the cumulative plastic strain variable p defined as 7.13 According to the state laws (7.10), we obtain 7.14 The regularization operator then follows from the combination of the previous constitutive equations with the balance equation (4.2): 7.15 The specific choice leads to a regularization operator that is linear with respect to Lagrangian coordinates: 7.16 which involves the Laplacian operator Δ0 in the isotropic case, i.e. , 7.17 The impact on hardening can be seen by choosing, as an example, α=p, according to (7.13), as an internal variable in (7.12). The dissipation potential Ω can be chosen in such a way that the residual dissipation takes the form 7.18 with , where f is the yield function. As a result, the yield stress R is given by the following enhanced hardening law: 7.19
(b) Eulerian formulation
In the Eulerian formulation, the free energy density is taken as a function of instead of , according to (7.2): , so that the state laws for generalized stresses become 7.20 The arguments of the free energy mix the invariant quantities and the observer-dependent variable . Galilean invariance then requires ψ to be isotropic with respect to .
The constitutive choice (7.12) is now replaced by 7.21 A quadratic potential ψ∇ is necessarily of the form , for objectivity reasons, so that the regularization operator involves the Laplace operator Δ w.r.t. Eulerian coordinates: 7.22 If the same viscoplastic yield function as in the previous subsection is adopted, the hardening rule is enhanced as follows: 7.23 This therefore yields a finite strain generalization of Aifantis strain gradient plasticity model .
(c) Formulation using the local intermediate configuration only
In the two previous formulations, Lagrangian or Eulerian generalized strain variables were mixed with the elastic strain variable attached to the intermediate local configuration, as the arguments of the free energy function. It is possible to assign the free energy function with a consistent set of arguments solely attached to the intermediate configuration. For that purpose, a generalized strain and generalized stresses are now defined on the intermediate local configuration: 7.24 and 7.25 The power density of internal forces expressed w.r.t. the intermediate local configuration takes then the form 7.26 where is still given by equation (7.7). To establish this expression, the following relation was used: 7.27 The dissipation rate density measured w.r.t. the intermediate local configuration is then 7.28 The free energy density function is chosen as . As such, it is invariant w.r.t. change of observer. The Clausius–Duhem inequality is now derived as 7.29 This expression reveals the existence of a generalized Mandel tensor, , conjugate to the plastic deformation rate, that is a function of the classical Mandel stress and of microdeformation related stress and strain. Positivity of dissipation is ensured by the choice of the following state laws and plastic flow and hardening rules: 7.30 and 7.31 provided that the dissipation potential displays suitable convexity properties with respect to both arguments.
The yield criterion is taken as a function where the ΠMeq is an equivalent stress measure based on the generalized Mandel stress tensor. Choosing α=p, where p is still given by equation (7.13), the residual dissipation takes the form 7.32 Note that the contribution in the generalized Mandel stress acts as a size-dependent kinematic hardening component which comes in addition to isotropic hardening represented by R. This is a specific feature of the model formulation w.r.t. the intermediate configuration.
As an example, a typical form of the free energy density function based on constitutive variables defined on the intermediate configuration, and hyperelastic laws are 7.33 and 7.34 These generalized stresses can be inserted into the balance equation 7.35 This provides the form of the regularization operator: 7.36 A quadratic dependence of the contribution leads to the following linear relationship between and . However, in that case, the regularization operator (7.36) is nonlinear and does not involve a Laplace operator, even in the isotropic case : 7.37 with and . As a result, the hyperelastic relationships for the higher order stresses are not linear w.r.t. to the associated strain gradient measures.
The construction of regularization operators presented in this work is based on the introduction of strain or damage-like micromorphic degrees of freedom in addition to the displacement vector and of their gradients into the Helmholtz free energy density function of the constitutive model. The combination of a new balance equation for generalized stresses and of the micromorphic constitutive equations generates the wanted regularization operator. Within the small strain framework, the choice of a quadratic potential w.r.t. the gradient term provides the widely used Helmholtz operator whose regularization properties are well demonstrated in the literature: smoothing of discontinuities at interfaces and boundary layers in hardening materials, and description of finite width localization bands in softening materials. The thermodynamic theory also predicts the form of the coupling between the standard and micromorphic variables. When the micromorphic degrees of freedom are related to total strain tensor or equivalent strain measures, the coupling arises in an extended elasticity law. For microplastic variables, the elasticity law is unaffected whereas the hardening law is modified.
Non-quadratic potentials w.r.t. the gradient term were proposed including power and logarithmic laws motivated by recent advances in strain gradient plasticity. They provide new strongly nonlinear differential operators whose mathematical properties were not explored and are largely unknown. They were formulated considering anisotropy, as required for applications to crystalline metals and damaging composites for instance.
The presented extensions to finite deformations show that the regularization operator cannot be postulated in an intuitive way. It is rather the result of a constitutive choice regarding the dependence of the free energy function on the gradient term. Purely Lagrangian and Eulerian formulations are straightforward and lead to Helmholtz-like operators w.r.t. Lagrangian of Eulerian coordinates. Two alternative standard procedures of extension of classical constitutive laws to large strains, widely used in commercial finite-element codes, have been combined with the micromorphic approach. In particular, the choice of local objective rotating frames leads to new nonlinear regularization operators that are not of the Helmholtz type. Three distinct operators were proposed within the context of the multiplicative decomposition of the deformation gradient. A new feature is that a free energy density function depending on variables solely defined with respect to the intermediate isoclinic configuration leads to the existence of additional kinematic hardening induced by the gradient of a scalar micromorphic degree of freedom.
Note that the results obtained for the micromorphic theory with additional degrees of freedom are also valid for gradient theories (gradient plasticity or gradient damage) once an internal constraint is imposed linking the additional degrees of freedom to strain or internal variables. This amounts to selecting high values of parameter Hχ or introducing corresponding Lagrange multipliers. The analysis was essentially limited to scalar micromorphic degrees of freedom for the sake of simplicity, even though tensorial examples were also given. Scalar plastic microstrain approaches suffer from limitations like indeterminacy of flow direction at cusps of the cumulative plastic strain in bending, for instance, see [60,97,98]. Those limitations can be removed by the use of tensorial micromorphic variable (microstrain or microdeformation tensors). The micromorphic approach is not limited to the gradient of strain-like, damage or phase field variables. It can also be applied to other internal variables, as demonstrated for hardening variables in [38,99].
It remains that the regularization properties of the derived nonlinear operators are essentially unknown, except through examples existing in the mentioned literature. For instance, the Eulerian and Lagrangian variants of the Helmholtz-type equation for scalar micromorphic strain variables have been assessed recently in  giving the advantage to the latter, based on finite element simulations of specific situations. The regularizing properties of more general operators should be explored in the future from the mathematical and computational perspectives in order to select the most relevant constitutive choices that may depend on the type of material classes.
It may be surprising that the constitutive theory underlying the construction of regularization operators for plasticity and damage, mainly relies on the enhancement of the free energy density function instead of the dissipative laws. It is in fact widely recognized that plastic strain gradients, e.g. associated with the multiplication of geometrically necessary dislocations, lead to energy storage that can be released by further deformation or heat treatments. However, the limitation to the enhancement of free energy potential is mainly due to the simplicity of the theoretical treatment and to the computational efficiency of the operators derived in that way. Dissipative higher order contributions remain to be explored from the viewpoint of regularization, as started in . Constitutive models of that kind are already available for plasticity, damage and fracture [101,102].
The author has no competing interests.
Part of this study was carried out in the framework of project MICROMORFING (ANR-14-CE07-0035-03) funded by the Agence Nationale de la Recherche (ANR, France). This support is gratefully acknowledged.
The author thanks his colleagues and PhD students at Centre des Matériaux Mines ParisTech for many fruitful discussions.
Appendix A. Notations for higher order gradients of a scalar field
Notations are introduced for the normalized gradient of a scalar field ϕ corresponding to a micromorphic variable, and also to a phase field: , where ℓc is a characteristic length introduced for normalization purposes. The direction of the gradient is given by the vector A 1 More generally, an (anisotropic) norm of the gradient is introduced, , involving a definite positive symmetric second-order tensor (dimensionless). The derivative of the anisotropic norm with respect to is denoted by and computed as A 2
The gradient of the anisotropic norm is calculated in a similar way: A 3 The derivative w.r.t. and the gradient of the anisotropic direction are computed as A 4 and A 5 In the isotropic case, , and .
- Received October 30, 2015.
- Accepted March 29, 2016.
- © 2016 The Author(s)
Published by the Royal Society. All rights reserved.