On configurational forces for gradient-enhanced inelasticity

In this paper we discuss how configurational forces can be computed in an efficient and robust manner when a constitutive continuum model of gradient-enhanced viscoplasticity is adopted, whereby a suitably tailored mixed variational formulation in terms of displacements and micro-stresses is used. It is demonstrated that such a formulation produces sufficient regularity to overcome numerical difficulties that are notorious for a local constitutive model. In particular, no nodal smoothing of the internal variable fields is required. Moreover, the pathological mesh sensitivity that has been reported in the literature for a standard local model is no longer present. Numerical results in terms of configurational forces are shown for (1) a smooth interface and (2) a discrete edge crack. The corresponding configurational forces are computed for different values of the intrinsic length parameter. It is concluded that the convergence of the computed configurational forces with mesh refinement depends strongly on this parameter value. Moreover, the convergence behavior for the limit situation of rate-independent plasticity is unaffected by the relaxation time parameter.


Introduction
Strain gradient effects in metals become particularly important in situations when the deformation localizes within a B Dimosthenis Floros florosd@chalmers.se 1 Department of Industrial and Materials Science, Chalmers University of Technology, 412 96 Gothenburg, Sweden narrow band, whose width is small as compared to a representative, problem-specific, length (e.g. the size of the plastic zone ahead of a crack-tip). That such localization zones can occur has been shown clearly in experiments, see e.g. the torsion tests on thin copper wires in Fleck et al. [1], and the indentation tests on single-and poly-crystal copper in Nix and Gao [2]. Although gradient effects are related to characteristics of the crystal lattice (microscopic scale), see Ashby [3], attempts in the literature during the last three decades focus a great deal on the derivation of constitutive theories that can be fully defined by macroscopic measures (such as the gradient of the plastic strain weighted by a material length, i.e. l s [ p ⊗ ∇]). Important contributions towards such phenomenological theories are found in Aifantis [4] and Zbib and Aifantis [5].
However, capturing the "right" physics related to size effects is not the sole reason for employing the gradient theory. In fact, regularization of otherwise non-smooth gradient fields in order to provide robust numerical solutions is another motivation behind the use of non-local theories in general, see Bažant and Jirásek [6], and more specifically gradient theories, see Svedberg and Runesson [7] and Liebe and Steinmann [8]. The phenomenological strain gradient theory introduced by Fleck and Hutchinson [9] is used in Xia and Hutchinson [10] in order to derive closed-form solutions for the near-tip fields, whereby it is required that the near-tip fields obtained for strain gradient plasticity should tend to the HRR fields from local theory (see Hutchinson [11]) at a distance from the crack-tip that is larger than the chosen internal length. In addition, Chen et al. [12], adopting the gradient-enhanced theory in Fleck and Hutchinson [13], compared the validity of the asymptotic near-tip fields with full-field numerical solutions.
Configurational (or material) forces are the energetic driving forces (sensitivities) pertinent to evolving configurations in mechanics. Examples are the thermodynamic force on a material inhomogeneity, cf. Eshelby [14], or the energy release rate in fracture mechanics, cf. Griffith [15]. Computation of such forces for inelastic materials most often involves the evaluation of spatial gradients of internal variables 1 . In turn, the way this evaluation takes place is determined by the "structure" of the variational problem at hand. Namely, in conventional displacement-based formulations, values of the internal variables are obtained from the solution of the time-discrete evolution equations at the integration points, i.e. only discrete values are available. In order to compute the spatial gradient in such a case, it is natural to utilize some sort of smoothing strategy. Such a strategy introduces further discretization errors, especially in crack problems involving steep gradients, see Tillberg et al. [16] and Özenç et al. [17]. In Özenç et al. [17] and Näser et al. [18], it is shown that pathindependence of the J -integral for inelasticity is achieved only by including the material dissipation inside the contour.
In Menzel et al. [19], it is shown that shifting from a conventional displacement-based variational formulation to a mixed displacement-internal variable formulation adds nothing to the accuracy of the measured configurational forces-related quantities, while the conventional formulation results in computationally more efficient schemes. Another type of mixed formulations is studied in Liebe et al. [20], where the displacements are coupled to an isotropic continuum damage variable.
In the present contribution, configurational forces are computed based on a gradient-enhanced constitutive theory via a mixed-dual variational formulation. A previously derived thermodynamically consistent definition of configurational forces, developed by Tillberg et al. [16] is thereby used. The coupled primary fields are the displacements along with a so-called micro-stress 2 field, the latter being the stress measure which is energy-conjugated to the spatial gradient of the internal variables. In this way, two interrelated goals are accomplished. Namely, smooth gradient fields are obtained, which together with the nodal values of the gradient field 3 , obtained from the solution of the primary problem at hand, provide a rational basis towards computation of configura- 1 The evaluation of spatial gradients of internal variables is performed in order to compute the so-called material dissipation part of the total configurational force, see e.g. Tillberg et al. [16]. 2 The respective field, which is properly introduced in Sect. 2, is known as "micro-stress" in the context of crystal plasticity, see Gurtin [21]. The relation between the micro-stress and a Peach-Koehler type of force (see Maugin [22]) is also shown in Gurtin [21]. The Peach-Koehler force may be viewed as the configurational force acting on an elastic dislocation. 3 The fact that nodal values of the gradient field are known already from the solution of the primary problem essentially means that no nodal smoothing of internal variables from computed values at the integration points is needed at the post-processing for the computation of configurational forces. tional forces that is more insensitive to the mesh resolution. The gradient effects in the developed model are scaled by an intrinsic length, which effectively acts as an internal regularization parameter.
Earlier work on configurational forces in the context of gradient-enhanced modeling (damage) concern general computational aspects, Frankenreiter et al. [23], as well as the particular issue of mesh-adaptive procedures, Welschinger and Miehe [24]. Theoretical aspects of configurational forces for non-local elasto-plasticity are also discussed in Stumpf et al. [25], with application to the evolution of fracture process regions.
The paper is organized as follows: in Sect. 2, the gradientenhanced mixed-dual variational formulation is presented. A brief outline of the configurational motion problem in the local and the gradient-enhanced constitutive settings is given in Sect. 3. Two numerical examples are investigated in Sects. 4 and 5, whereby the solution of the primary problem is validated, and the characteristics of configurational forces are examined. Further reflections on results from the primary problem and the configurational motion problem are discussed in Sect. 6. Finally, conclusions are presented in Sect. 7.

Primal constitutive setting
As a prototype model, we consider perfect viscoplasticity of the Bingham type which is enhanced by an appropriately formulated gradient term. Assuming that the gradient effect is energetically decoupled from the local effect, we propose the free energy density where [u] = [u ⊗ ∇] sym is the strain 4 , p is the plastic strain (the single internal variable) and g[ p ] = p ⊗ ∇ is the gradient variable. Henceforth, we denote p as "internal variable", although it is a field variable in the non-local setting. Square brackets [ ] denote operational dependence of the argument. Furthermore, we consider linear elasticity and linear gradient hardening where the modulus tensor E e = 2G I sym dev + K I ⊗ I represents isotropic elasticity, whereas H g l 2 s may be viewed as a gradient hardening modulus. In this paper, H g is taken as constant, while the internal length l s may be viewed as a regularization parameter, i.e. l s = 0 corresponds to the classical local model.
The "energetic" stresses that are conjugated to , p and g are given as To complete the model, we introduce the (dual) dissipation potential φ * in terms of the "dissipative" stress κ di representing Bingham viscoplasticity, as follows: where η(F) is the overstress function and F is the (quasistatic) yield function. ( ) dev denotes the deviatoric part, whereas = 1 2 ( + | |) defines the Macaulay bracket. The material parameters are t * (relaxation time), σ y (yield stress) and σ c (characteristic stress). The flow rule, i.e. the evolution equation for p , is thus given aṡ

Balance equations in primal form
The strong format of the physical problem under quasi-static conditions and in the absence of body forces reads: find the fields u(x, t), p (x, t) and κ di (x, t) that satisfy Equation (7) is the standard equilibrium under quasi-static conditions, while Eq. (8) pertains to the "micro-force" balance equation, see e.g. Gurtin [21]. Fig. 1 Body occupying the domain , surface tractions t p , microtractions p p and normal n. Dual partitioning of the boundary is introduced as We may reformulate Eq. (8) as κ di = σ + ξ · ∇ = "σ red ", whereby it appears that −ξ · ∇ plays the role of backstress in kinematic hardening for purely local theory.

Balance equations in mixed-dual form
Employing a Legendre transformation, we may replace g by ξ as an independent variable in the expression of the free energy. The pertinent transformation reads where " . . ." denotes the scalar product between third order tensors.
A semi-dual free energy denoted ϕ( , p , ξ ), can now be defined as For the specific gradient-enhanced viscoplastic material model of Bingham type introduced in Eqs. (1) and (2), the semi-dual free energy ϕ( , p , ξ ) in Eq. (17) becomes Moreover, the constitutive equation for g(ξ ) is obtained from Eqs. (16) and (18) as As a result of this change of variables, we replace the set of Eqs. (7)-(9) by the following mixed-dual format: Find the fields u( The time-discrete version of Eqs. (20)-(23) is obtained using the Backward-Euler time integration rule in the time interval (t n , t n+1 = t n + t]. Thereby, all the field variables are considered known at time t n , while their values at time t n+1 are sought.

Mixed-dual weak form
The time-discrete variational form of Eqs. (20)-(23) is given next. Firstly, the appropriate trial and test spaces are introduced, where L 2 ( ) is the Hilbert space of square integrable functions, H 1 ( ) and H 1 div ( ) are vector spaces with gradients and divergence, respectively, in L 2 ( ). The mixed-dual weak form of the time-discrete problem can now be stated as follows: Find u ∈ U, p ∈ K, ξ ∈ X, and κ di ∈ K that solve where Here, the upper left index n indicates a quantity at t = t n , whereas the lack of temporal index refers to the sought value at t = t n+1 . The boundary terms on the RHS of Eqs. (29) and (31) read The appropriate boundary conditions, as introduced already in Eqs. (10)-(13), become apparent from the structure of the linear forms Eqs. (33) and (34). Alternatively, the resulting weak form Eqs. (29)-(32) may also be derived as the stationarity conditions of the incremental pseudo-potential, which can be further simplified by elimination of κ di from Eq. (37). Thus, we obtain the "local" residual where l (u) (δu) and l (ξ ) (δξ ) are defined in Eqs. (b) Compute R u (u (k) , ξ (k) ; δu) and R ξ (u (k) , ξ (k) ; δξ ).

Configurational motion
The continuum mechanics setting applied in the present work follows Runesson et al. [26]. The standard quasi-static equilibrium problem can be stated as that of analyzing the motion of a body initially occupying the (material/Lagrangian) domain . In response to the acting physical forces, the body will deform into the spatial (Eulerian) domain ω(t), which is time-dependent in general. The mapping between Lagrangian and Eulerian domains is y = ϕ(x, t), while the pertinent deformation gradient is defined as F def = ϕ ⊗ ∇. The small strain setting used in Sect. 2 can be identified upon defining u = y − x, whereby = [F − I] sym .
Configurational motion can now be introduced as a change in the material domain with time, i.e. = (t). To this end, a reference absolute/fixed configuration ξ , parametrized in ξ 5 , is introduced, from which it is possible to describe changes in the initial domain (e.g. a crack advance), see Fig. 2. Finally, a composite mapping may be defined between the fixed and Eulerian configurations, y = x + u =φ(ξ , t).
We are now in the position to describe the configurational motion as t x, often denoted configurational velocity, where t = ∂ ∂t | ξ is the time derivative 6 for fixed coordinates ξ ∈ ξ . Finally, in order to define discrete forces in subsequent sections, we may parametrize the configurational velocity in terms of a discrete velocity D t a and a scalar field W (x) as t x = W (x)D t a. Here, we evaluate W (x) = 1 at the crack-tip, defining the crack propagation. The invariance to how W is chosen inside will be discussed in Sect. 3.2.

Rate of mechanical dissipation and configurational forces
For the isothermal case, the total mechanical dissipation for a given finite body that is subjected to loading such that it deforms and, at the same time, undergoes configurational changes may be expressed as We shall henceforth assume that the configurational velocity is non-zero only on the part of ∂ where the tractions vanish. Now, adopting the free energy of a gradient-enhanced dissipative material in the small strain setting according to Eq. (1), thereby introducing the functional dependence ψ = ψ( , p , g), we may split the total mechanical dissipation in Eq. (46) in two parts, i.e.
In Eq. (47), the Eshelby energy momentum tensor is identi- Following the commonly adopted model that internal variables are "frozen" in the material domain, i.e. p (x, t) and g(x, t) are unaffected by the configurational velocity t x, cf. Runesson et al. [26], we obtain the sensitivities for a change in the configurational velocity δ t x. The directional derivative of D with respect to an arbitrary change δ t x, followed by the parametrization δ t x = W (x)δD t a, result in the generalized configurational force The last equality in Eq. (51) follows from the constitutive relations in Eqs. (3) and (4) and the definition of In fracture mechanics, the crack driving force acting on a crack-tip can be obtained by projecting G onto the tangential direction of the undeformed crack. The pertinent projection is commonly known as energy release rate (see Griffith [15]), and is equivalent to the J-integral introduced by Rice [27] (for linear elasticity). In contrast to the latter scalar-valued quantities, the benefits of formulating the configurational motion problem in terms of a vectorial quantity (such as the configurational forces) are summarized in Steinmann [28].
Remark 2 (Invariance of W ) Any configurational motion that only lives in the interior of the domain (i.e. vanishes on ∂ ) will not affect the dissipation. This can also be identified by the balance of material/configurational forces, since As a result, the choice of W (x) in the interior of will affect only the numerical approximation of G, and not its continuous counterpart. We note that the decomposed contributions G CONF and G MAT may still vary with respect to the choice of W (x) in the interior of , however, their sum is path-independent.

Computation of configurational forces
Consider a body with an embedded discontinuity (e.g. a crack) occupying the domain , as shown in Fig. 3. We will examine the computation of the material dissipation part of the configurational force, G MAT , for local theory and a displacement-based variational formulation, as compared to the gradient-enhanced mixed-dual formulation proposed here. The computation of the configurational part, G CONF , from Eq. (50) is performed in the same way for the considered variational formulations.
Assuming local theory and a displacement-based variational formulation, and employing the definition of configurational forces in Eq. (50), we conclude that b mat takes the form where the fields σ and p are computed at the integration points. Thereby, in order to compute g[ p ] at the integration points, a scheme of the form is necessary, see Fig. 3. In Eqs. (54) and (55), N a are nodal shape functions and m is the total number of nodes in a finite element discretization. For the general inelastic case, the equations for the variation of p in the near-tip region are not known. Conventionally, it is assumed that it suffices to obtain the nodal values p a via an L 2 -projection of the known values of p at the integration points onto the space of first or higher order polynomials. In the case of inelastic response, nodal smoothing techniques such as the aforementioned L 2projection may turn out insufficient for capturing the steep gradients occurring in the near-tip region, contributing to the "pathological" FE-mesh sensitivity reported in e.g. Tillberg et al. [16].
In the context of the gradient-enhanced mixed-dual variational formulation described in Sect. 2, the vectorial quantity b mat in Eq. (51) takes the form where the last equality follows from the constitutive relations in Eqs. (3) and (19) and the assumption that the material properties are uniform within the body. In contrast to Eq. (53), a pronounced advantage of the mixed-dual formulation proposed here is that all quantities required for the computation of Eq. (56) are known already from the solution of the primary problem for (u, ξ ). Namely, [u], ξ and ξ ⊗ ∇ are known explicitly from the FE-approximation, whereas p is evaluated at the Gauss-points. In the following sections, the computation of configurational forces for the gradient-enhanced formulation is evaluated in two problems: (a) virtual expansion of a smooth interface, see Sect. 4, and (b) virtual extension of a discrete crack, see Sect. 5. The main goal of these investigations is to examine whether the proposed gradient-enhanced mixed-dual formulation provides sufficient regularity for the computation of configurational forces. In addition, we seek to identify conditions that allow for the configurational forces obtained from the proposed mixed-dual formulation to approach those obtained by a "conventional" displacementbased variational formulation based on local theory.

Problem description
The problem of a square plate with a centric hole is examined, see Fig. 4a. The plate is subjected to prescribed vertical displacement of the upper edge u 2 = 0.1 mm, while the displacement is fully restrained at the bottom edge, i.e. u 1 = u 2 = 0.0 mm. The prescribed displacement varies linearly in time from 0 mm to its peak value from t = 0 s to t = 0.05 s, which is the final time of the simulation. Plane strain conditions are assumed. The continuous domain is discretized in triangular elements with piecewise linear approximations of displacements u and micro-stresses ξ , see Appendix C. The FE-mesh is structured, graded towards the hole, as shown in Fig. 4b, so that the steep gradients closest to the hole region are adequately resolved. The material parameters that enter the gradient-enhanced perfect viscoplasticity model of the Bingham type in both examples outlined in Sects. 4 and 5 are given in Table 1.
In accordance with the discussion on boundary conditions occurring from Eq. (8) in Sect. 2.2, two possible choices are examined. Namely, either "hard" boundary conditions are  imposed on the plastic strain, p = 0, or, "free" boundary conditions are prescribed on the micro-traction, p = ξ · n = 0, along the entire boundary. In the first case, the pertinent external source term in Eq. (34) vanishes. In the second case, only the components of ξ parallel to the outward pointing unit normal n on the boundary vanish, for the boundaries that are parallel to one of the axes of the coordinate system introduced in Fig. 4a. As for the ξ -degrees-of-freedom that "live" at the hole boundary, imposing "free" boundary conditions necessitates the appropriate transformation of the micro-stress tensor ξ from the global coordinate system of Fig. 4a to a local system at each point on the hole boundary, with the "x-axis" being the outward pointing unit normal n in the specific point.

Model validation and influence of boundary conditions
In Fig. 5, the distribution of the independent 7 plastic strain components over the entire domain is illustrated for a sufficiently fine mesh, at the last incremental loading step (i.e. for u 2 = 0.1 mm, see Fig. 4a), for the "hard" boundary con-ditions, p = 0 and for finite value of the internal length, l s = 0.02 m. From the pertinent images, it is evident that the plastic strain components tend to vanish along the entire boundary of the domain, as it is expected from the specific choice of weakly prescribed boundary conditions on the plastic strain.
A mesh convergence study of the plastic strain components In particular, the prescribed "hard" boundary conditions, p = 0, are being better captured as the mesh resolution becomes finer. It should be noted that since the plastic strain field is known only at the integration points, the curves shown in the latter figure correspond to smoothed nodal values obtained via L 2 -projections of the known values of the plastic strain at the integration points. Details on the computation of the L 2 -projections can be found in e.g. Hinton and Campbell [29]. The mesh convergence of the plastic strain field along the same line for a model based on local (visco)plasticity theory is also included for comparison, see Fig. 7. It is verified that, unlike the case of the gradient-enhanced model under "hard" boundary conditions, the plastic strain field converges to a finite non-zero value in the case of local (visco)plasticity theory. This is also known from theory, since the hole boundary comprises a smooth interface, thus no singularity effects are present. Such singularity effects occur in cases where sharp cracks are embedded in an elastic (see Irwin [30]) or elastoplastic (see Rice and Rosengren [31]) continuum.
At this point, it should be noted that the curves labeled as "local" throughout this study have been obtained via a local viscoplastic material model (i.e. the free energy of the model is ψ = ψ loc ( , p ), where ψ loc ( , p ) is defined  1)). The use of a local viscoplastic model imposes certain characteristics and limitations, which are discussed further in Sect. 6. The most prominent is that for such a formulation, no boundary condition on the plastic strain or the micro-traction can be, nor need to be, imposed.
The influence of the boundary conditions, pertinent to micro-force balance equation (8), is examined next. The resulting plastic strain fields for the specific "hard" and "free" boundary conditions employed in this example are depicted More specifically, from the normalized plastic strain fields in Fig. 8, it becomes evident that, due to the "hard" boundary conditions, p = 0, the respective curves for the gradientenhanced models are "weakly" forced to vanish close to the hole boundary. However, for the model that follows the local theory, the pertinent curve converges to a finite non-zero value, as shown in Figs. 7 and 8. In sharp contrast, the curve for the lowest value of the length l s = 0.0002 m in Fig. 9, which corresponds to "free" boundary conditions, p = 0, approaches to a large extent the behavior of the "local" curve. Thus, it may be argued that "free" boundary conditions on the micro-traction p is a more appropriate choice, if the goal is to approximate the response of the local theory.

Analysis of convergence with respect to the FE-mesh size
In this section, we investigate the mesh sensitivity of the configurational forces. The configurational forces are computed according to Sect. 3, based on the "equilibrium fields" defined by the solution of the primary problem described in Sect. 4.1. The material model defined in Sect. 2 and the respective material parameters outlined in Table 1 are also adopted for these numerical investigations.
The particular configurational motion considered here is uniform expansion of the hole, as shown in Fig. 10. The position vectors that describe the relevant motion are depicted in Fig. 10a, and the length of those vectors is illustrated via isolines in Fig. 10b. In turn, the isolines correspond to values of the scaling function W in Eq. (50) at the nodes that define the isolines. The total energy release rate for a (virtual) increase of the hole radius can thus be obtained from a radial projection of Eq. (50) integrated along the circumference. More specifically, where e r is a unit vector defining the radial direction from the center of the hole. In Figs. 11 and 12, the energy release rates are depicted for "hard" boundary conditions, p = 0, and "free" boundary conditions, p = 0, respectively. We note that G is the total energy release rate (as defined in Eq. (57)), which is comprised of the configurational part, G CONF , and the material dissipation part, G MAT . The mesh sensitivity of the global response in terms of the total vertical reaction force at the upper edge due to the prescribed displacement of that edge, see Fig. 4a, for different values of the internal length, is depicted in Fig. 11f for "hard"  1) and (2). As regards energy release rates, the gradient-enhanced model proposed here results in computable quantities for both the total rate, G, as well as the configurational and the material dissipation parts, G CONF , and G MAT , see Fig. 11ac, respectively. Moreover, we introduce the decomposition G MAT = G MAT,term1 + G MAT,term2 , where the two terms stem from the two terms of b mat in Eq. (56). The mesh sensitivity of these terms can be viewed in Fig. 11d, e. That the magnitude of these terms is almost identical is noteworthy.
In addition, the curves shown in Fig. 11a-e do not appear to be ordered with respect to the value of the internal length, neither for the total, nor for the configurational and the material dissipation parts. The interpretation of this state of affairs is addressed later in this section. The sensitivity of the total vertical reaction force, due to prescribed displacement of the upper edge, is studied in Fig. 12f, for "free" boundary conditions, p = 0. As was also observed in Fig. 11f, the global response of the gradient-enhanced model converges to the one from local theory as the internal length l s decreases towards zero.
Focusing on the mesh sensitivity of the energy release rates for "free" boundary conditions, p = 0, we note that convergence is obtained for all the considered quantities, namely G, G CONF , and G MAT , see Fig. 12a-e. As is also the case for "hard" boundary conditions, p = 0, convergence appears to be "slower" for smaller values of the internal length. This behavior is expected, since the models with lower values of the internal length result in less regularized (smooth) gradient fields, p and g, which, in turn, requires finer mesh resolution for resolving these gradients.
Concerning the response for "free" boundary conditions, p = 0, it is clear from Fig. 12b, c that the depicted curves are ordered with respect to values of the internal length. This is in accordance with the discussion on the choice of boundary conditions for the gradient problem in Sect. 4.2. Namely, in Fig. 8 it was shown that prescribing "hard" boundary conditions, p = 0, results in the solution of a different primary problem in terms of e.g. the plastic strain field, as compared to the one for local theory. In contrast, prescribing "free" boundary conditions, p = 0, results in plastic strain fields that are able to capture the response from local theory for decreasing l s , see Fig. 9. Moreover, according to the theory and as verified by the local model employed in this study, the plastic strain at the hole boundary converges to a finite non-zero value, see Fig. 7. The facts that (a) only by prescribing "free" boundary conditions, it is possible to capture the plastic strain distribution from the local theory near the hole boundary, and (b) the maximum value of the configurational motion takes place at the hole boundary (see Fig. 10), provide sufficient explanation as to the particular ordering of the curves in Fig. 12b, c.

Analysis of convergence with respect to the internal length
In Sect. 4.3.1, the mesh sensitivity of the energy release rates was examined, where it was shown that the proposed gradient-enhanced mixed-dual variational formulation leads to computable rates. It was also shown that finer mesh resolution is required for the models with lower values of the internal length to converge. The convergence error of G MAT with respect to the ratio of the internal length to the characteristic element length, l s / h, is depicted in logarithmic scale for both axes, in Fig. 13a, b for "hard" boundary conditions, p = 0, and "free" boundary conditions, p = 0, respectively. In the latter figures, the relative error in G MAT is computed according to the expression No of displacement dof where the "exact" solution, denoted here as G MAT h , is taken as the value of G MAT for the finest mesh, see Figs. 11c and 12c. The observed rates are approximately quartic for "hard" boundary conditions, p = 0, and quadratic for "free" boundary conditions, p = 0.
The variation of G, as the internal length l s tends to zero, is illustrated in Fig. 14a, b for "hard" boundary conditions, p = 0, and "free" boundary conditions, p = 0, respectively. The depicted curves correspond to values of G for the finest mesh, see Figs. 11a and 12a for l s = 0.0002 m, 0.002 m, 0.02 m. The value of G from local theory is also marked with dashed line. It is observed that, for "free" boundary conditions, the behavior of the local theory is approached by the gradientenhanced model, as l s tends to zero. In the case of "hard" boundary conditions, the limit value of G from the gradientenhanced model converges to a finite value, which is much larger than the value of G from local theory. No of displacement dof 5 Numerical example: discrete crack

Problem description
The problem of a single edge-cracked specimen under plane strain loading is examined next. The specimen is subjected to a prescribed displacement u 2 = 2.4 mm of the upper edge, while both degrees-of-freedom (dof) at the bottom boundary are fixed, i.e. u 1 = u 2 = 0.0 mm, see Fig. 15a. The simplest FE-approximation for both fields (u, ξ ) is employed here, i.e. triangular elements with piecewise lin-ear approximations of displacements u and micro-stresses ξ , see Appendix C, Fig. 22. Five unstructured FE-meshes are used in this investigation, graded towards the cracktip (see Fig. 15b), featuring a characteristic element length of 34.0 mm, 19.0 mm, 10.5 mm, 6.3 mm, 4.7 mm. Regarding boundary conditions pertinent to micro-force balance equation (8), the same alternatives are examined here, as for the problem with the smooth interface outlined in Sect. 4. Namely, either "hard" or "free" boundary conditions are imposed along the whole boundary.

Analysis of convergence with respect to the FE-mesh size
Based on the "equilibrium fields" obtained by the primary problem described in Sect. 5.1, in principle, the same quan- tion, as depicted in Fig. 16b. These projections correspond to the J-integral for elasticity, see also Sect. 3.2. Stated in terms of configurational motion, unit crack advance that is tangential to the crack-tip direction is examined. The pertinent configurational motion is scaled here, linearly, according to Fig. 16a. The values depicted in the figure denote the values of the scalar-valued function W in Eq. (50). For the straight edge-crack studied here, the projection of the configurational force in the direction tangential to the undeformed crack-tip corresponds to the tangential component of this force, as measured in the coordinate system introduced in Fig. 15a. The perpendicular component of the configurational force is not examined here, as it is shown in the literature to be path-dependent, even under elastic material response, see e.g. Brouzoulis and Ekh [32]. The mesh sensitivity of the energy release rate G and its constituents, G CONF and G MAT , is examined in Figs. 17 and 18 for "hard" and "free" boundary conditions, respectively. The global response of the gradient-enhanced model for different values of the internal length is compared with the response obtained from local constitutive theory in Fig. 17f, for "hard" boundary conditions, p = 0, and Fig. 18f, for "free" boundary conditions, p = 0. In both cases, the global Footnote 9 continued coordinates. In turn, in the present context, the material domain may vary if crack propagation is applicable. In such a case, the undeformed crack-tip direction would be continuously updated along with the Lagrangian domain. response of the gradient-enhanced model approaches that of local theory as the internal length is reduced.
It can be argued from Figs. 17a and 18a, that the proposed gradient-enhanced mixed-dual variational formulation results in a computable quantity G . Convergence for finite values of the length l s is obtained for a relatively coarse FE-discretization, while smaller values require a finer mesh resolution. This behavior is similar to the mesh sensitivity of the energy release rates for the smooth interface problem outlined in Sect. 4.3.
Examining the constituents of G , namely G CONF and G MAT , we observe that the pertinent curves are ordered with respect to the values of the internal length, see Figs. 17b, c and 18b, c for "hard" boundary conditions, p = 0, and "free" boundary conditions, p = 0, respectively. For the depicted values of l s , the behavior of the gradient-enhanced model approaches that of local theory for decreasing values of the internal length, irrespective of the choice of "hard" or "free" boundary conditions. However, the employed mesh resolutions do not allow for drawing a definite conclusion regarding the value of G for infinitesimal l s , see Sect. 5.2.2.
Unlike the example of the smooth interface (cf. Fig. 11d,  e), the constituents of G MAT , i.e. G MAT,term1 and G MAT,term2 , are not of the same order of magnitude, see Figs. 17d, e and 18d, e.

Analysis of convergence with respect to the internal length
How the relative error of G MAT reduces with the ratio of the internal length to the characteristic element length, l s / h, is studied in Fig. 19a, b, for "hard" boundary conditions, p = 0, and for "free" boundary conditions, p = 0, respectively. The relative error is computed according to Eq. (58). The value of G MAT for the finest mesh is taken as the "exact" value G MAT ,h in this analysis, see Figs. 17c and 18c. It is observed that the error reduces with quadratic or higher rate with l s / h, irrespective of the choice of boundary conditions. This is similar to the rate of convergence of G MAT for the smooth interface problem for "free" boundary conditions, p = 0, cf. Fig. 13b.
The variation of G , as l s tends to zero, is illustrated in Fig. 20a, b, for "hard" boundary conditions, p = 0, and for "free" boundary conditions, p = 0, respectively. The value of G from local theory is also marked with dashed line. The depicted curves pertain to values of G for the finest mesh, see also Figs. 17a and 18a for l s = 0.02 m, 0.2 m, 1.0 m. It is therefore emphasized that, for infinitesimal l s and especially for the local theory, the depicted values of G in Fig. 20 have not fully converged for the considered FE-mesh sizes in Figs. 17a (and 18a). From Fig. 20a, b, it is observed that the choice of "hard" or "free" boundary conditions has little influence on the limit value of G obtained from the gradientenhanced model.

Rate-independent plasticity as a limit case of viscoplasticity
In this section, the convergence properties of the configurational forces, as computed by the proposed gradientenhanced mixed-dual formulation, are investigated for decreasing values of the relaxation time, t * = 0.1 s, 0.056 s, 0.01 s, 0.001 s. The pertinent decrease corresponds to a transition of the gradient-enhanced prototype model of perfect viscoplasticity defined in Sect. 2.1, from a viscous to a rate-independent (t * → 0) material response. All the other material parameters given in Table 1 are kept constant for this numerical investigation. In Fig. 21a, b, G CONF and G MAT are depicted for fixed value of the internal length, l s = 0.02 m, and for "hard" boundary conditions, p = 0. Similar convergence behavior is observed, irrespective of the choice of the relaxation time. In addition, the behavior of the gradient-enhanced model approaches that of rate-independent plastic response, as the relaxation time diminishes. As expected, lower values of t * lead to higher amount of material dissipation, represented by G MAT , due to a more pronounced inelastic response.

Discussion
In the present contribution, it has been shown that it is possible to obtain gradient fields with sufficient regularity for the computation of configurational forces, see e.g. Figs. 11, 12, 17 and 18. As expected, this comes at the "cost" of requiring a fine mesh resolution for small value of the internal length. It was shown in Sect. 5.2.3 that the proposed formulation results in similar convergence behavior for fixed values of the internal length l s irrespective of the value of the relaxation time t * . Next, the possible choices regarding the value of the internal length are reviewed. In particular, if the internal length is viewed as a model parameter that needs to be calibrated, its value is determined by the characteristic dimension of an underlying substructure (such as the size of the grains in the near-tip region). If, on the other hand, the internal length is viewed as a regularization parameter, the limit case l s → 0 is of interest. In such a case, an optimal value of the internal length is chosen, which comprises a compromize between (a) the proximity to the gradient fields obtained from the pure local constitutive theory and (b) the regularity of the gradient fields at hand.
The behavior of G as l s tends to zero for the problem of the smooth interface is depicted in Figs. 14a and 14b, for "hard" boundary conditions, p = 0, and "free" boundary conditions, p = 0, respectively. The limit value of G is different, depending on the choice of boundary conditions related to micro-force balance equation (8), while considering that no such choice can be made for the "standard" displacementbased formulation. More specifically, the limit value of G from the gradient-enhanced model converges to the value obtained from local theory, only for "free" boundary conditions. For "hard" boundary conditions, G converges to a finite value which is much larger than the value from local theory. The highlighted difference in the limit value of G for the two choices of boundary conditions is linked to the effect this choice has on the resulting plastic strain fields already at the primary problem, cf. Figs. 8 and 9 for "hard" and "free" boundary conditions, respectively. In addition, the choice of "hard" boundary conditions seems as the least physically sound (for both the considered problems), since under that choice, the hole (or the crack-tip) is forced to remain elastic, when it is known that it is there where most of the stress concentration takes place.
The behavior of G as l s tends to zero for the single edge-cracked specimen is depicted in Fig. 20a, b, for "hard" boundary conditions, p = 0, and "free" boundary conditions, p = 0, respectively. Unlike the problem of the smooth interface, the choice of "hard" or "free" boundary conditions appears to have little influence on the limit value of G . The fact that the choice of boundary conditions does not affect the limit value of G is attributed to the extent of the plastic zone size in the specific problem, which corresponds to an infinitesimal region closest to the crack-tip. The corresponding region in the problem of the smooth interface spans a much larger portion of the domain, near the hole boundary.
Aspects of computational efficiency are discussed next. Specifically, choosing a sufficiently large value for the internal length ensures convergence even for coarse FE-meshes, see e.g. Fig. 18a. The latter choice of length implies that gradient effects are influencing a much larger portion of the domain. This results in numerical models which are computationally more "expensive" as compared to formulations with less pronounced gradient effects for the same FE-discretization. Therefore, the gains originating from the coarser FE-meshes required for finite internal lengths to be resolved, should be weighted by the loss stemming from the "heavier" primary problems at hand for finite values of the internal length.
For the proposed mixed formulation, additional dof per node for the gradient field are introduced in addition to the displacement dof. Hence, to reduce the computational cost, more flexible schemes may be constructed. For example, as regards single edge-cracked specimen, it may be sufficient to introduce the proposed gradient-enhanced formulation only within a region closest to the crack-tip, where the gradients are "steeper". In the remaining part of the domain, "conventional" displacement-based variational formulation and local constitutive theory may be used. Regarding the boundary conditions on the interface between the gradient-enhanced and local formulations, the different extremes "free" or "hard" could be considered. The specific boundary conditions may serve as an indication of the distance from the crack-tip that the gradient-enhanced formulation needs to be considered. Beyond that distance, the influence of the boundary conditions should clearly vanish.
Furthermore, the proposed mixed-dual formulation is comprised of displacements and micro-stresses, which may have different orders of magnitude. This implies that the global tangent stiffness matrix may become ill-conditioned. Proper scaling of the stiffness matrix should be a good remedy, however, further detailing is outside the scope of this paper.

Conclusions
Configurational forces have been computed via a gradientenhanced mixed-dual variational formulation. An internal length was used as a regularization parameter. A significant advantage is that no nodal smoothing of internal variables is required at the postprocessing. Investigation of the mesh sensitivity related to the computed configurational forces has highlighted that the pertinent formulation results in computable quantities for both types of problems examined in the paper: a smooth interface and an edge crack. The discretization error was found to decrease with quadratic or higher rate with respect to the ratio of the internal length to the characteristic element length l s / h.
The configurational forces computed by the proposed gradient-enhanced formulation were also compared to configurational forces obtained from a "conventional" displacement-based variational formulation for local constitutive theory. It was found that, (a) for the smooth interface, the proposed gradient-enhanced model resembles the behavior of the local theory as the internal length tends to zero, only when homogeneous boundary conditions are imposed on the micro-traction, and (b) for the edge crack, the choice of homogeneous boundary conditions on the plastic strain or the micro-traction has little influence on the response of the gradient-enhanced model at the limit of vanishing internal length.
In conclusion, adopting the gradient-enhanced model allows for numerically well-behaved evaluation of the configurational force, at the cost of solving a more "expensive" problem than that for classical local theory. Depending on boundary conditions, the converged result may (or may not) be strongly affected by the introduced material length scale.
where the terms ∂ p /∂ and ∂ p /∂χ are determined from the linearization of the local residual, see Eqs. (71) As is evident from Eq. (38), the plastic strain p is implicit function of and χ. Therefore, the linearization of the local residual R L with respect to perturbations d , dχ and d p provides explicit expressions for ∂ p /∂ and ∂ p /∂χ which enter the expressions for the linearized residuals Eqs. (59)-(62). The pertinent linearization takes the form Performing the differentiations described in Eq. (66) gives where C, A and B are determined by Considering Eq. (67), we are now able to state explicitly the expressions that determine ∂ p /∂ and ∂ p /∂χ , provided that p is symmetric and the "operator" 11 C is invertible. 10 Clearly, Eq. (63) is exactly the algorithmic tangent stiffness in case of local plasticity theory, defined by ξ = 0. 11 The term "operator" is used here to highlight the fact that the fourth order tensor C can be viewed as an operator that transforms d p into the equivalent expression on the RHS of Eq. (67).

Reduced-Voigt format of mixed-dual problem in continuous and discrete spatial domains
In Sect. 2, a gradient-enhanced mixed-dual variational formulation was proposed, in terms of displacements and micro-stresses, (u, ξ ), in small strain setting. The resulting weak form is defined by Eqs. (39), (40). If standard Voigt notation is adopted, the latter equations (neglecting the external work) may be rewritten as where denotes the Voigt format of tensor , i.e. a column matrix containing all cartesian components of the tensor. In Voigt notation, the following relations between the tensors and the respective conversion matrices may be written: where A σ , A , A p , A χ , A g , A ξ (A ξ = A g ) are conversion matrices that convert the stored values of σ v , v , p,v , χ v , g v and ξ v , from a "reduced-Voigt format", to their full (3D) tensor equivalents in Voigt notation, which would read 9 × 1 and 27×1 components for a second-and a third-order tensor, respectively. Explicit forms of all the entries in Eqs.
where δ v ≡ v [δu] and δχ v ≡ χ v [δξ ]. Without loss of generality, the reduced-Voigt format for a plane strain problem is discussed next. In view of this assumption, the non-zero independent components in σ and reduce to 4 and 3, respectively. As regards p , it should be noted that the pertinent tensor is symmetric and deviatoric, therefore the independent components that need to be stored in each integration point reduce to 3. The latter properties of p simplify also the computation of g (and consequently also of ξ ), a fact that stems from the definition of g, namely g = p ⊗∇, which yields g (and ξ ) symmetric with respect to its first two "legs". Therefore, it suffices to store only 6 components of g in each integration point (cf. 27 components for the full (3D) tensor and 10 for a plane strain problem without accounting for the facts that p is symmetric and deviatoric). Symmetry of p yields χ symmetric, thus 4 non-zero components for χ need to be accounted for, under the considered plane strain setting. In turn, the independent components of σ , , p , χ, g, ξ are stored in vectors σ v , v , p,v , χ v , g v and ξ v , respectively. The latter vectors with stored values at the integration points enter Eqs.
g v = [ g 111 g 221 g 121 g 112 g 222 g 122 ] T , where γ 12 = 12 + 21 , while noting that 12 = 21 . Combining Eqs. (82)-(87) and (90)-(95) fully defines the entries in the conversion matrices A σ , A , A p , A χ , A g , namely and Next, the discretization of the field variables u and ξ is performed, along with the equations that fully define the continuous system, i.e. Eqs. (88)-(89). The primary variables are approximated according to where N u (x), N ξ (x) are matrices of nodal shape functions, and a u , a ξ are nodal values. Note that under the proposed gradient-enhanced mixed-dual formulation, FEapproximations of different order for u(x) and ξ (x) may generally be used. Similarly, the test functions δu, δξ are approximated as which is in accordance with Galerkin's method, in which the same shape functions are employed for both the primary variable and the test function. In addition, the fields (x) and χ (x) are discretized according to For the plane strain problem discussed here, the simplest possible finite element for the considered mixed FE-problem is the 3-node triangular element, i.e. with linear approxi- ]. The number of independent nodal ξ -dof that need to be accounted for is 6, in accordance with the preceding discussion pertaining to the independent components of g.
The global residual equations in reduced Voigt format, Eqs. (88) and (89), may also be written in discrete form by replacing the continuous fields u, ξ v , δu, δξ v , v and χ v from Eqs. (98)-(103). The pertinent equations read where in the frames of the plane strain setting and the considered number of independent components in the relevant tensors, the FE-matrices B u , B ξ and N ξ take the form with n being the total number of nodes.